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On the modeling of car passenger ferryship 
design parameters with respect to selected 
sea-keeping qualities and additional 
resistance in waves 


Tomasz Cepowski, Assoc. Prof. 
Szczecin Maritime University 


ABSTRACT 


y 
e second part of the investigations design solutions were searched for by applying the single- 
and multi-criterial optimization methods. The multi-criterial optimization was performed by 
using Pareto method. Such approach made it possible to present solutions in such form as to allow decision 
makers (shipowner, designer) to select solutions the most favourable in each individual case. 


This paper presents the modeling of car passenger ferryship design parameters with respect 
to such design criteria as selected sea-keeping qualities and additional resistance in waves. 
In the first part of the investigations approximations of selected statistical parameters of 
design criteria of ferryship were elaborated with respect to ship design parameters. The 
approximation functions were obtained with the use of artificial neural networks. In the 


Keywords: sea-keeping qualities; roll-on - roll-off ferryship; rolling; motion sickness index; lateral accelerations, additional 
resistance in waves; ship design parameters; modeling; artificial neural networks; optimization; Pareto method 


INTRODUCTION 


In ship design process design solutions which fulfil both 
economic criteria and technical limitations are searched for. 
Economic criteria are consisted ofa set ofrequirements imposed 
by shipowner, among which a suitable internal capacity and 
service speed of the ship can be numbered, that significantly 
influences operational profitability of the ship sailing on a given 
shipping line. Obtaining the assumed service speed by the 
ship depends a.o. on parameters and operational conditions 
of propulsion system as well as value of total hull resistance 
to motion. Since total resistance of ship hull consists of a.o. 
additional resistance in waves associated with ship sailing in 
heavy weather conditions, the additional resistance can be used 
as a design criterion for certain types of ships. 

And, technical limitations are consisted of a.o. appropriate 
stability, unsinkability or hull structural strength. For certain 
types of ships insensibility to weather conditions, i.e. the so- 
called good sea-keeping qualities are more and more often 
used as a design criterion. In such case among sea-keeping 
qualities wave effects can be numbered in the form of rolling 
and secondary wave effects (resulting from waving and rolling) 
to which a.o. accelerations, slamming, shipping of water on the 
deck and sickness motion index, belong. 

The presented investigations were aimed at determination 
of optimum values of ship design parameters with respect 
to assumed economic criteria and technical limitations. 
The investigations were performed on the example of 


preliminary design of a car passenger ferry ship. The type of 
ship is characteristic of specific operational and design problems 
such as internal capacity, stability (unsinkability), speed or 
insensibility to weather conditions. Therefore for the ship in 
question selected sea-keeping qualities and additional resistance 
in waves at a given ship capacity were assumed to be design 
criteria. In the first phase of the investigations approximation 
functions of the above mentioned criteria were determined on 
the basis of main geometrical parameters of ship hull. In the 
second phase values of the design parameters were determined 
for which the assumed design task were characterized by the 
best merits as regards the assumed criteria. 


METHOD 


The following design task was formulated in the 
investigations: to find the vector of independent variables, 
X=X (X X, --- X,), which minimizes or maximizes functions 
of partial targets under limitations given for car passenger ferry 
ship sailing in heavy weather conditions. 

The following data were assumed: 

— independent variables: L/B, B/d, CB, CWL, where: L, B, d 
- ship length, breadth and draught, CB - block coefficient 
of hull underwater part, CWL - waterplane coefficient; 

— limitations: 

e theoretical volumetric displacement V = 17500 m° 

e L/B ratio within the range of 5.17 + 7.64 

e B/d ratio within the range of 3.22 + 4.46 
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e CB within the range of 0.6 + 0.64 

e CWL within the range of 0.8 + 0.85 

e initial lateral metacentric height GM = 1 m 

e ship speed v = 5 kn; 

— design criteria: 

e the motion sickness index MSI (acc. ISO 2631/3) [9] 
for the wave encounter angle B = 120° (in the reference 
system of: 180° — head wave, 0° — following wave) 

e the additional resistance in waves R for B = 180° 

e the roll angle o for B = 30° 

¢ the lateral acceleration on the car deck, at, acc. [11] for 
B= 30°. 


In the investigations was assumed a conventional wave 
spectrum consisted of wave energy spectral density functions 
acc. ITTC for waves of the significant heights H = 6 m and 
characteristic periods T in the range from 6 to 14 s, (Fig.1). 
This made it possible to eliminate impact of the characteristic 
period on ship responses and to decrease number of independent 
variables in the design task. 

10 


[= 


Conventional wave energy 


spectral density function [m?/s] 


03 04 05 06 0.7 08 09 1 11 12 13 
Wave frequency [1/s] 
Fig. 1. Conventional wave energy spectral density function, 
ITTC spectrum, H = 6 m, T = 6 + 14 s 


In the first phase of the investigations functions of partial 
targets were determined depending on independent variables. 


APPROXIMATION OF DEPENDENT 
VARIABLES 


Functions of partial targets, which are necessary to solve 
the design task, are to be expressed in the form of analytical 
functions approximating the assumed dependent variables 
(i.e.: motion sickness index, additional resistance in waves, roll 
angles and lateral accelerations). To perform approximations 
the method presented in Fig. 2 was assumed. 


Calculation of model values of dependent variables: 
y= Y: (Xb X2 Xn )i=1..4 


Calculation of the functions f approximating the 
dependent variables y: 
f, 
Y; Xp X2 Xn) f (Xp, Xq,---Xn) 


Fig. 2. Algorithm of approximation of the dependent variables y, 
where: X „ X, ... X,— independent variables (design parameters), 
y, — model values of dependent variables, 

J, — searched for approximating functions, i=1..4. 
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To determine the functions f, in the equation (1) it is 
possible to make use of statistical methods or those based on 
artificial neural networks. The application of artificial neural 
networks to approximation of unknown relations belongs to 
mathematical numerical methods dealing with the so-called 
„artificial intelligence” and finds wider and wider use in 
shipbuilding [1, 2, 5, 6]. 

To calculate model values the set of 24 design variants 
of car passenger ferry ship of different hull forms shown in 
Tab. 1, was assumed. Model values were presented in the form 
of irregular wave statistical values of the wave energy spectral 
density function shown in Fig. 1. 

Tab. 1. Model set of independent and dependent variables appearing in the 
above assumed parameters of ship motions and waving, where: CB — block 
coefficient of hull underwater part, CWL — waterplane coefficient, L — ship 
length, B — ship breadth, d — ship draught, MSI — motion sickness index, 
a,— lateral acceleration on car deck, @— significant amplitude of ship roll 
angles, R — additional resistance in waves 


Independent variables Dependent variables 


97.20 | 10.82 
88.40 | 9.56 
93.90 | 7.98 
91.50 | 9.61 
92.30 | 11.60 
91.40 | 7.50 
88.20 | 9.62 
66.80 | 9.57 
81.30 | 6.32 
79.50 | 6.27 
80.00 | 8.69 
79.30 | 7.66 
78.20 | 8.73 
92.80 | 10.51 
85.80 | 7.54 
80.80 | 7.84 
78.70 | 10.01 
75.20 | 6.80 
87.30 | 7.31 
79.60 | 9.07 
79.00 | 6.30 
71.20 | 5.54 
90.70 | 7.96 


878.50 
1150.60 
1164.00 
856.60 
796.40 
985.90 
712.00 
893.00 
933.60 
1068.10 
843.80 
970.30 
1108.90 
852.50 
1116.00 
956.30 
879.10 
1106.20 
878.90 
882.10 
1090.00 
798.90 
595.80 
775.50 


z 
& 
i 
R 
> 
í= 
g 
3 
a 
1 
2 
3 
4 
5 
6 
7 
8 
9 
10 
11 


— = =... m 
NYDN BW WN 


The values of dependent variables given in Tab. 1 were 
calculated by means of exact numerical methods with the use 
of SEAWAY software which is based on planar flow theory. 
Exactness tests of the software, presented in [4], indicate 
a high accuracy of calculations. For calculations of additional 
resistance in waves the Gerritsma-Beukelman method was 
used [3]. 

Approximation functions of the dependent variables: 
MSI, ọ, R and a, were determined by using artificial neural 
networks and presented analytically in the form of Eqs. (1), 
(2), (3) and (4): 


R VAW 


2 


| -0.423 -0.886 1.422 0.046 -0.087 -1.134 -0.560 


l 
m CWL,LB,Bid]xS+P)xA-B) © C- a =% 
MSI = 
2 
where: 
MSI — motion sickness index [%] 
L — ship length 
B — ship breadth 
d — ship draught 


CB — block coefficient of hull underwater part 
CWL -— waterplane coefficient 
A — matrix of weight values: 


-0.581 1.019 -2.326 -0.015 1.460 -0.044 
0.379 -0.268 -0.469 -0.021 0.310 -1.453 
-0.924 0.480 -1.928 -1.033 1.071 0.487 
-0.321 0.842 1.699 -0.686 -0.422 0.652 
-0.283 -0.254 2.564 1.254 1.030 0.340 


S — matrix of coefficients: 


0 0 
0 2083 0 0 0 
0 0 2041 0 0 
0 0 0 064 0 
| 0 0 0 0 O81 | 


— vector of threshold values: [-1.583 -0.131 -1.882 0.366 0.508 -2.022] 

— column vector of weight values: [1.15 -0.71 -2.80 1.59 -1.88 1.75] 

— vector of displacement values: [-13.31 -16.75 -15.02 -3.29 -2.60] 
a,— coefficients of the following values: a, = -1.62, a, = -2.74, a, = 0.04 


l 
f + e (ICWL.CBICWL.LB.Bid]x$+P)xA-B) xC- Oy |— Oy 


a2 
where 
R — additional resistance in waves [kN] 
A — matrix of weight values: 


0.726 2.897 -2.504 0.570 1.932 -0.365 -0.164 
3.661 -4.110 0.244 0.431 3.790 0.275 -0.070 
-1.219 -0.517 1.854 -0.481 0.820 -0.656 -0.492 


S — matrix of coefficients: 


20.83 0 0 0 
0 20.41 0 0 
0 0 0.64 0 
0 0 0 0.81 


vector of threshold values: [2.220 -2.108 1.688 -0.967 -0.410 -1.842 0.633] 
column vector of weight values: [2.665 2.431 1.836 0.373 -2.001 -2.197 -0.728] 
vector of displacement values: [-16.75 -15.02 -3.29 -2.60] 


— coefficients of the following values: a, = -0.939, a, = -1.575, a, = 0.0022 


j| 
ers CWL.LB.Bid]xS+P)xA-B) * C- Op |— Oy 
a, = SSS 
t 


9) 
where: 
a, — maximum lateral acceleration on the car deck [m/s’] 
A — matrix of weight values: 
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(1) 


(2) 


(3) 


-0.874 0.408 0.046 0.313 0.001 0.863 -0.096 0.041 0.368 0.265 0.732 -0.753 -0.471 
1.062 -0.014 -0.680 0.739 0.903 -0.981 -0.583 -0.089 -0.835 0.239 0.035 -0.380 -0.550 
0.955 -0.358 -0.392 0.412 -0.270 0.704 -0.467 0.149 -0.867 0.060 0.801 -0.202 0.736 
-0.755 0.616 -0.266 0.527 0.320 0.507 -0.968 -0.821 0.870 0.378 0.261 -0.480 -0.089 
| 0.810 -0.399 -0.822 0.233 0.480 0.903 -0.424 0.241 -0.655 0.385 0.301 -0.656 0.420 | 
S — matrix of coefficients: 
| 2222 0 0 0 0 | 
0 21.47 0 0 0 
0 0 20.41 0 0 
0 0 0 0.64 0 
| 0 0 0 0 0.81 


— vector of threshold values: 
[0.368 -0.043 -0.045 -0.231 -0.109 0.821 -0.289 0.924 0.776 0.683 1.025 -0.397 -0.121] 
— column vector of weight values: 
[-0.923 0.021 -0.562 0.617 0.171 -0.052 0.356 -0.773 0.798 0.786 0.628 1.003 -0.989] 
P — vector of displacement values: [-13.31 -17.48 -15.02 -3.29 -2.60] 
dp 0,, a, — coefficients of the following values: a, = 0.193, a, = -2.18, a, = 1.67 


(([CB,CWL,CB/CWL,B/d]xS+P)xA+6.96)+1.049 


Q Ww 


OP cl? 


po (4) 
0.189 
where: 
b — significant amplitude of roll angles [°] 
A — column vector of weight values: [ 38.14 -30.37 -33.83 -0.54] 
S — matrix of coefficients: 
22.22 0 0 0 
0 20.83 0 0 
0 0 20.41 0 
0 0 0 0.81 _J 


P — vector of displacement values: [-13.31 -16.75 -15.02 -2.60] 


In Tab. 2 the above presented networks are described 
together with their selected statistical parameters. From the 
above given sheets it results that the functions (1), (2), (3) and 
(4) are characterized by a simple structure and relatively high 
accuracy. 


Tab. 2. Type, structure and selected statistical parameters 
of the elaborated neural networks 


Root mean square error, 
RMS, of: 


Variable 
network 


learning 


testing 


Type of network 
Correlation 
coefficient R 


Structure of neural 


MLP 
MLP 
MLP 


5x6x1 

4x7x1 

5x13x1 
4x1 


linear 


SELECTION OF OPTIMUM DESIGN 
PARAMETERS OF SHIP HULL 


e maximum of MSI 

e minimum of the additional resistance in waves, R 

e minimum of the acceleration on the car deck, a, 

* minimum of the roll angle 9. 

The functions (1), (2), (3), (4) were used to solve the problem. 
To find a design solution which satisfies the assumed criteria the 
single- and multi-criterial optimization methods were applied. 


SINGLE-CRITERIAL OPTIMIZATION 


In this part of the investigations minimum values of the 
functions (1), (2), (3) and (4) were determined. In the first 
phase the GRG2 (Generalized Reduced Gradient) program of 
linear optimization was used [7]. The results are presented in 
Tab. 3 and Fig. 3 and 4. 

The main disadvantage of optimization methods of both 
kinds (single- and multi-criterial) is that they provide only 
information on location of an optimum point. They do not give 
information on shape of objective function and it is not certain 
whether the achieved optimum solution is global [3]. Hence 
in the next phase of the investigations a set of all possible 
design solutions was determined and the results were analyzed 
against the assumed criteria. To the calculations the following 
increments of the design parameters were assumed: 


e AL/B = 0.01 m 

In this part of the experiment hull form parameters of car * AB/d = 0.01 m 
passenger ferry ship were searched for with respect to the * ACB = 0.01 
assumed criteria, i.¢.: e ACWL = 0.01. 
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Tab. 3. Optimum design solutions with respect to the assumed criteria 


min. of a, | min. of ọ | min. of R | min. of MSI | Criterion 


500 600 700 800 900 
R [kN] 
Fig. 3. Comparison of the design solutions with respect to the additional 
resistance in waves, R, and the motion sickness index MSI 


1000 1100 1200 


— 
in 


a, [m/s?] 
A 


5 6 7 8 9 


4 
> [deg] 


Fig. 4. Comparison of the design solutions with respect 
to the roll angle and the lateral acceleration a, 


The analysis results are illustrated in Fig. 5 through 11. In 
Fig. 5, 7, 9, 11 are presented the sets of minimum values of 
the variables (design criteria) dependent on a given criterion. 
From the sets merits of the design solutions result with respect 
to the assumed criteria. For instance from Fig. 5 it results that 
the design solutions characterized by small values of MSI index 
are also characteristic of large values of additional resistance 
and small values of roll angles and lateral accelerations. And, 
the design solutions characterized by large values of MSI index 
are also characteristic of a.o. small values of mean additional 


resistance in waves. In Fig. 6, 8, 10, 12 values of design 
parameters are presented depending on selected design criteria. 
On the basis of both groups of the diagrams the car passenger 
ferryship design parameters which satisfy decision maker with 
respect to assumed criteria, can be determined. 


su] % *[ 


č 


= added resistance 


roll 
80 g 


50 55 60 65 70 75 


(OLX) [ 


Fig. 5. Minimum values of dependent variables 
within the range of MSI = 53 + 88 % 


8 0.64 
7 
0.635 
6 small R 
THA large a, 
v’ large 0.63 Q 
a 4] largeR = 
4 small a, ie 
m3 small > 0.625 
ar l 
0.62 
1 
0 0.615 
50 55 60 65 70 75 80 85 90 
MSI [-] 
Fig. 6. Values of design prameters in function of MSI index, CWL = 0.80 
18 
16 
= 
14— 
= 
12 98, 
10 £ 


(or [su 


4 

2 

0 

500 600 700 800 900 1000 1100 
R [kN] 


Fig. 7. Minimum values of dependent variables 
within the range of R = 530 + 1060 kN 


large MSI 
large 
large a, 


min. MSI 
min. 
small a 


L/B [-]. Bid [-] 


0. 
500 600 700 800 900 
R [KN] 
Fig. 8. Values of design parameters in function 
of the additional resistance in waves R, CWL = 0.80 


1000 
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c 


GOTS) [s/t] *% [%4] 


$ [deg] 


Fig. 9. Minimum values of dependent variables 
within the range of ġ = 0.26 + 6.93° 


7 

6 

5 
a 
v 
= 4 
mM 
Gra mean MSI 
—3 small R 
a small a, 
2 


3 4 
[deg] 
Fig. 10. Values of design parameters in function of the roll angle ọ 


MSI [%]. x [deg] (x10) 


1.15 Lz 1.25 1.3 1.35 1.4 1.45 1.5 
a, [m/s*] 


Fig. 11. Minimum values of dependent variables 
within the range of a, = 1.19 + 1.46 m/s? 


small R 
large MSI 
large 


large R 
small MSI 
small > 


1.15 1.2 1.25 1.3 1.35 1.4 1.45 LS 
a, [m/s?] 


Fig. 12. Values of design parameters in function 
of the lateral acceleration a, CWL = 0.80 
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MULTI-CRITERIAL OPTIMIZATION 


The next step in the investigations was an attempt of using 
the mult-criterial optimization methods in order to determine 
an optimum design solution of ro-ro ferry ship. In such case the 
notion of optimum solution is not simple as in the case when 
one does not accept to compare different criteria to each other 
and then must propose such definition of optimality which 
respects integrity of each of the criteria [5]. In order to select 
an optimum solution in the analyzed case when sets of partial 
solutions do not have any common part, it is necessary to find 
an acceptable compromise between partial targets. Importance 
of particular partial targets can be expressed by means of the 
weight coefficient [3] and hence obtained design solutions 
depend on values of the weights. In such case a verdict is 
automatically given and relations between partial solutions are 
unknown. This is disadvantageous as decision maker has no 
possibility to assess all variants (or only a part of them) before 
a verdict is given. 

In the case of the ship hull parametric modeling for the 
offer design purposes it is sometimes important to know 
influence of partial criteria on design solutions. Then in the 
offer design its merits can be proposed from the point of view 
of different criteria and decision maker (shipowner) is able 
to choose the best variant from his point of view. Therefore 
in the assumed design task the notion of optimality in the 
sense of Pareto was used. The Pareto optimality concept 
does not provide any indications as to making choice of 
a final solution (out of P-optimum ones) . Decision on 
choice of optimum solution has to be taken by decision 
maker himself [5]. 

The optimality condition in the Pareto sense can be 
presented as follows: the vector x will be consider (partly) 
smaller than the vector y, if and only if the relation [5] is 
satisfied: 


(Vi)(x; < y;,)A Gi); < y;) (5) 


Then it is also said that the vector y is dominated by the 
vector x. If a given vector is not dominated by no other one 
then it is called the non-dominated vector. By analyzing Fig. 3 
one can state that the solutions A and B are not dominated and 
the solutions C and D are dominated by the solutions A and B 
(with respect to MSI index and additional resistance R). And, 
from Fig. 4 results that the soultions A and B are dominated 
by the solution D with respect to the criteria of the roll angle 
ọ and lateral acceleration a,. From the fact it results that in the 
assumed design task no single optimum solution with respect to 
all the criteria can be obtained but only the set of four solutions 
which are not mutually dominated, and any decision on choice 
of a solution is compromise. 

In Fig. 13 is presented the set of values of the dependent 
variables of all the optimum solutions - in the sense of Pareto 
- (i.e. those non-dominated) of the considered design task. The 
set is presented in function of the MSI index. As it turned out, in 
the obtained design solution the dependent variables achieved 
the following values: 

e the motion sickness index, MSI = 55 + 90 % 

e the additional resistance in waves, R = 614 +1164 kN 

e the significant amplitude of roll angles = 0.26 + 8.91° 
e the lateral accelerations a, = 1.22 + 1.91 m/s’. 


From the above given set optimum solutions can be selected 
with respect to arbitrary assumed criteria. In Tab. 4, 5 and 6 are 
presented example design solutions with respect to: 


— all the criteria which achieve maximum mean values, 
i.e.: MSI < 85.5%, R < 889 KN, p < 4.59°, a, < 1.57 m/s? 
(Tab. 4) 

— MsSland R values smaller than mean ones, as well as @ anda, 
values greater than mean ones, i.e.: MSI < 78%, R< 730 kN, 
 < 7°, a, < 1.7 m/s? (Tab. 5) 

— anda, values smaller than mean ones and MSI and R values 
greater than mean ones, i.e.: MSI < 90%, R < 1165 KN, 
 < 2°, a, < 1.39 m/s? (Tab. 6). 


== 
| 


Fig. 13. The set of all optimum design solutions in the sense of Pareto, 
MSI = 55 + 90 %, R = 614 +1164 kN, f = 0.26 + 8.91°, a, = 1.22 + 1.91 m/s? 


Tab. 4. Optimum design solutions in the sense of Pareto, 
MSI < 85.5%, R < 889 KN, ġ < 4.59°, a, < 1.57 m/s? 


3.57 |85.46) 883.65 | 2.97 | 1.59 


Tab. 5. Optimum design solutions in the sense of Pareto, 
MSI < 78%, R < 730KN, ġ < 7° a, 1.7 m/s? 


Tab. 6. Optimum design solutions in the sense of Pareto, 
MSI < 90%, R < 1165 KN, $< 2°, a, < 1.39 m/s? 


89.67 | 1164.2 


CONCLUSIONS 


— Both economic criteria and technical limitations constitute 
important factors in the process of single- and multi-criterial 
design optimization. To perform such optimization ship 
designers should be equipped with suitable tools in the 
form of calculation methods. Determination of additional 
hull resistance in waves or ship sea-keeping qualities by 
means of the presently used design algorithms is excessively 
complicated as to make it possible to perform optimization 
successfully. This results from various causes but the most 
important is that simple and simultaneously exact relations 
between hull form parameters and ship behaviour in waves 
are lacking. To predict ship responses in waves is only 
possible by using exact numerical methods for a given ship 
hull form. 


— A method which makes it possible to solve the problem is 
presented in this work. To approximate additional resistance 
in waves and selected sea-keeping qualities, artificial neural 
networks were used. Despite relatively low number of 
model variants and wide range of ship hull parameters the 
obtained approximations are characterized by high accuracy 
and simple structure. 


— The presented approximations may be applied to design 
analyses (e.g. such as in [2]) or as objective functions in 
various optimization methods of ship design parameters. 
In this paper an example of single- and multi- criterial 
optimization of car passenger ferry ship design parameters 
is presented. In the first phase ranges of optimum design 
parameters with respect to the assumed criteria taken into 
account separately, were determined. In the second phase 
sets of design solutions optimum in the sense of Pareto, 
were determined. In each of the phases the solutions were 
presented in such form as to give a decision maker (ship 
operator, designer) possibility to select variants most 
favourable to him. 
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ABSTRACT 


The article presents the results of the research project concerning the process of formation of the tip vortices 
shed from hydrofoils of different geometry in different flow conditions. Three hydrofoils resembling the 
contemporary marine propeller blades have been selected for the study. The experimental part of the project 
consisted of the LDA measurements of the velocity field in three cross-sections of the vortex generated by 
the hydrofoils in the cavitation tunnel. The numerical part of the project consisted of calculations of the 
corresponding velocity field by means of three computer codes and several selected turbulence models. The 
comparative analysis of the experimental and numerical results, leading to the assessment of the accuracy 
of the numerical methods, is included. 


Keywords: marine propulsors; vortex generation; numerical methods; experimental techniques 


INTRODUCTION 


One of the most important problems in the advanced design 
of marine propulsors is the possibly accurate determination of 
the concentrated vortex structures generated by the lifting foils 
forming the propulsor. These vortex structures usually lead to 
formation of the cavitating kernels, which are responsible for 
the cavitation erosion of the propulsor and the rudder, for the 
generation of pressure pulses leading to vibration of the ship 
structure and for the emission of intensive hydro-acoustic 
signals. Development of the methods for the accurate prediction 
of the above described phenomena at the design stage would 
lead to the complete elimination of their negative consequences 
or at least to limiting them to the acceptable levels. 

The rapidly developing methods of the Computational 
Fluid Mechanics are generally regarded as unable to predict 
accurately the intensity of the concentrated vortex structures 
generated by lifting foils, hence they are also unable to predict 
correctly the conditions for cavitation inception in these vortex 
structures. It is believed that the structures of the discrete 
representation of the fluid domain typically employed in these 
methods are unsuited for the cases with high velocity secondary 
transverse flows and that the standard turbulence models are 
developed and calibrated mostly for typical flat plate boundary 
layers, differing significantly from the vortex dominated flows. 
The objective of the research project presented in this paper 
is to develop the rules for generation of the discrete meshes 
especially suited for vortex dominated flows and to select the 
turbulence models especially effective in such flows. 


This objective is achieved by the combined experimental 
and numerical research. Three variants of the hydrofoil were 
selected for the project. Their outline and section profiles 
resemble those of the typical contemporary marine propeller 
blades. The radial distribution of pitch typical for a propeller 
blade was reduced to zero, thus producing the foil with standard 
distribution of the hydrodynamic loading along the span, 
denoted N for neutral. As the distribution of the hydrodynamic 
loading along the span of the foil may influence significantly 
the process of formation of the tip vortices, two more foils 
were designed: one with positive linear twist along the span, 
denoted L for loaded, and another with negative twist along 
the span, denoted U for unloaded. The absolute value of the 
twist angle in both cases was equal to 2.5 degrees at the tip, 
corresponding to the typical values for marine propeller blades. 
The details of the experimental and numerical parts of the 
project are described below, with the following comparative 
analysis of the experimental and numerical results. 


EXPERIMENTAL MEASUREMENTS 
OF THE VELOCITY FIELD AROUND 
THE CONCENTRATED VORTICES 


The comprehensive descriptions of the experiments with 
three hydrofoil models are presented in [1, 2, 3]. All three 
hydrofoil models were manufactured in aluminum alloy using 
the numerically controlled milling machine, which ensured 
high accuracy (Fig. 1). The span of the hydrofoil models 
was selected equal to 400 mm, which located the trailing tip 
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vortex approximately in the centre of the measuring section 
of the cavitation tunnel of CTO S.A. The hydrofoils were 
installed in the measuring section of the cavitation tunnel. 
Using the dynamometer, which enabled both measurements 
of the lift and drag forces and the accurate setting of the angle 
of attack (Fig. 2). The nominal velocity of flow in the tunnel 
was set to 4.0 [m/s], which ensured the turbulent flow over 
the hydrofoil. 


Fig. 1. The model of the hydrofoil N manufactured in aluminum alloy using 
the numerically controlled milling machine 


Fig. 2. The model of the hydrofoil fixed in the measuring section 
of the cavitation tunnel at CTO S. A. 


rrr 


Fig. 3. The laser Doppler anemometer (LDA) ready for measurements 
of the velocity field around the vortex behind the hydrofoil 


Each of the three hydrofoils was in turn set at three angles of 
attack: +2.5; 0.0 and —2.5 degrees. These values correspond to 
the typical range of variation of the angle of attack encountered 
on marine propellers. The measurements of the velocity field 
behind the hydrofoils were conducted by the laser Doppler 
anemometer (LDA — Fig. 3). The LDA system was capable of 
measuring two velocity component simultaneously: the axial 
component u and the transverse (vertical) component v. The 
measurements were conducted in three surfaces perpendicular 
to the vortex axis, located at 10 mm, 70 mm and 300 mm behind 
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the hydrofoil (Fig. 5). Apart from the velocity field, the variation 
of the lift and drag forces on the hydrofoils as the functions of 
the angle of attack were also measured. 


NUMERICAL PREDICTION 
OF THE VELOCITY FIELD AROUND 
THE CONCENTRATED VORTICES 


The comprehensive details of the numerical prediction of 
the velocity field in the vortex wake behind the hydrofoils are 
described in [1, 2, 4, 5]. Three numerical codes were used: the 
commercial programs Fluent and Comet and an indigenous 
program Solaga developed at CTO S.A. All three programs 
are based on the Reynolds Averaged Navier Stokes Equations, 
converted into algebraic equations by means of the Finite 
Volume Method. The computational domain had the same 
cross-section as the cavitation tunnel, i.e. 800 x 800 mm and it 
extended 1500 mm in front of the hydrofoil axis and 3000 mm 
behind the hydrofoil axis. This domain was divided into the 
block-structured grid having approximately 1.250.000 finite 
volumes — see Fig. 1. The size of the smallest finite volumes 
close to the hydrofoil surface was approximately 0.02 mm. The 
following turbulence models were used in calculations: 

e In the program Fluent 6.3, using the scheme MUSCL 
(Monotone Upstream-Centered Schemes for Conservation 
Laws) for the convection terms of the transport 
equations: 

- Spalart-Allmaras 

- k-epsilon RNG 

- k-omega SST 

- Reynolds Stress Model (RSM) 

e Inthe program Comet: 

- k-epsilon 

- k-omega 

- k-omega SST 
e Inthe program Solaga — Spalart-Allmaras 


Fig. 4. The discrete mesh for numerical calculations on the hydrofoils 
surface (top) and on the external surfaces of the flow domain (bottom) 


Fig. 5. The location of the three control surfaces both for measurements 
and calculations of the velocity field behind the hydrofoil 


The calculations reproduced the experimental conditions as COMPARISON OF THE EXPERIMENTAL 


closely as possible and they concentrated on determination of AND NUMERICAL RESULTS 
the axial and transverse (vertical) velocity components in the 
three cross-section of the vortex wake, respectively 10 mm, The results selected for comparison refer to two cross- 
70 mm and 300 mm behind the hydrofoil trailing edge, as sections of the vortex wake, namely 10 mm and 70 mm behind 
shown in Fig. 5. the hydrofoil. At both sections, the results are presented for 

0.37 

0.38 

0.39 
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0.41 

N 0.42 

0.43 
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0.47 
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0.04 0.02 0 -0.02 0.04 0.02 0 -0.02 
Y Y 
Comet (k-epsilon) Comet (k-omega SST) 


Fig. 6. The comparison of the measured (top) and calculated (below) axial component of the velocity 
at the control surface 10 mm behind the hydrofoil L at the angle of attack equal +2.5 degrees 
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these turbulence models which demonstrated the best ability to 
reproduce the experimental results. All presented results were 
obtained for the hydrofoil L set at the angle of attack equal to 
+2.5 degrees, i.e. for the case of the highest hydrodynamic 
loading at the tip of the hydrofoil. Consequently, the most 
intensive tip vortex was expected in this condition. Fig. 6 
shows the measured and calculated axial velocity component 

0.37 


0.38 
0.39 
0.40 
0.41 
0.42 
0.43 
0.44 
0.45 
0.46 
0.47 


at the cross-section nearest to the hydrofoil trailing edge. This 

velocity component is the reflection of the viscous wake behind 

the hydrofoil, generated in the boundary layer on its surface. 

The analysis of the velocity distributions shown in this figure 

leads to the following conclusions: 

— calculations using all programs and turbulence models give 
the qualitatively correct picture of the viscous wake behind 


0.04 0.03 0.02 0.01 0.00 -0.01 -0.02-0.03-0.04 


SRL TET aD 


Y 
Comet (k-epsilon) 


0.02 


Y 
Comet (k-omega SST) 


Fig. 7. The comparison of the measured (top) and calculated (below) axial component of the velocity 
at the control surface 70 mm behind the hydrofoil L at the angle of attack equal to +2.5 degrees 
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the hydrofoil, which agrees quite well with the results of — this over-prediction of the minimum velocity in the wake 
is much stronger for the RSM and k-epsilon turbulence 


the measurements, 
the calculated boundaries of the viscous wake are more 


models than for the k-omega SST turbulence model, 


produced very similar results using the same turbulence 


sharply defined than those determined experimentally, — itis interesting to notice that both Fluent and Comet have 
the calculated minimum velocity in the viscous wake 
is markedly smaller than that detected in the course of model k-omega SST. 
experiments, 
037 
0.38 
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0.39 
0.50 
0.40 
0.41 0.00 
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Fig. 8. The comparison of the measured (top) and calculated (below) transverse component of the velocity 
at the control surface 10 mm behind the hydrofoil L at the angle of attack equal to +2.5 degrees 
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Fig. 7 shows the analogical comparison of the measured 
and calculated axial component of the velocity at the cross- 
section located at 70 mm behind the hydrofoil trailing 
edge. As could be expected, at this location the viscous 
wake behind the hydrofoil was already significantly 
dispersed and damped due to the action of turbulence 
and fluid viscosity. This effect was correctly qualitatively 


0.02 0 -0.02 -0.04 


Y 
Comet (k-epsilon) 


reproduced in all calculations. However, again the k-omega 
SST turbulence model gave the best agreement with the 
results of experimental measurements, both in terms of the 
spatial distribution of the wake and the minimum velocity. 
Similarly as for the first cross-section, the k-epsilon and 
RSM turbulence models again over-predicted the minimum 
velocity in the wake. 


1.00 
0.50 
0.00 
-0.50 
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Fig. 9. The comparison of the measured (top) and calculated (below) transverse component of the velocity 
at the control surface 70 mm behind the hydrofoil L at the angle of attack equal to +2.5 degrees 
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In Fig. 8 the comparison of the calculated and measured 
transverse (vertical) component of the velocity in the wake 
behind the hydrofoil is presented for the cross-section located 
10 mm behind the hydrofoil trailing edge. This velocity 
component may be regarded as the measure of the intensity of 
the tip vortex. The analysis of the velocity distributions leads 
to the following conclusions: 

— calculations using all programs and turbulence models 
give the qualitatively correct picture of the intensity of the 
tip vortex, which agrees quite well with the results of the 
measurements, 

— even the small details of the spatial distribution of the 
vertical velocity component are very well reproduced in all 
calculations (cf. the shape of the red field in the pictures), 

— again the k-omega SST turbulence model seems to give the 
results closest to those obtained experimentally. 


Fig. 9 presents the distribution of the horizontal component 
of the velocity at the cross-section 70 mm behind the hydrofoil. 
Contrary to the picture shown in Figs. 6 and 7, now there is 
almost no turbulent and viscous dispersion effect of the vortex 
visible between Figs. 8 and 9. 

Simply the vorticity is still being convected towards the 
centre of the tip vortex. Consequently, the distributions shown 
in Figs. 8 and 9 are very similar to each other. This phenomenon, 
detected experimentally, is confirmed by all calculations. In 
this case no significant differences between the respective 
turbulence models may be noticed. 


CONCLUSIONS 


After a thorough analysis of the presented results the 
following general conclusions may be drawn: 

— the state of the art Computational Fluid Dynamics methods 
are fully capable of predicting the existence, location, 
spatial extent and intensity of the tip vortices generated by 
hydrofoils, 

— moreover, they are also capable of predicting correctly the 
intensity and consequences of two mutually counteracting 
processes of dissipation and concentration of vorticity along 
the vortex, 

— there are visible differences between the performance 
of different turbulence models and selection of the 
appropriate model is an important matter — in case of vortex 
dominated flows the k-omega SST model seems to the most 
promising, 


— there is definitely a room for improvement of the accuracy 
of numerical prediction of the tip vortex flows through 
modification of the discrete grid structure in the flow 
domain, especially in the direct vicinity of the vortex. 


The further stages of the research project, concerning 
the improvement of the performance of the CFD methods in 
prediction of the vortex dominated flows will be the subject 
of a separate article. 
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ABSTRACT 


The present study focuses on the simulation of two dimensional unsteady cavitating flows. For simulation 
of unsteady behaviors of cavitation which have practical applications, the development of unsteady PISO 
algorithm based on the non-conservative approach is utilized. For multi-phase simulation, single-fluid 
Navier-Stokes equations, along with the volume fraction transport equation, are employed. The bubble 
dynamics model is utilized to simulate phase change between vapor and liquid phases of the water. Unsteady 
simulation of cavitation around NACA66(MOD) and supercavitation around a flat plate normal to flow 
direction are performed to clarify accuracy of presented model. Numerical results and comparisons with 
experimental data are provided. The accuracy is good, and it is possible to apply this method to more 
complex shapes. 


Keywords: numerical simulation; cavitation; psio algorithm; flat plate; unsteady flow 


INTRODUCTION 


When local static pressure of liquid falls below the 
corresponding saturated pressure, the phase of fluid changes 
from liquid into vapor. This phenomenon is named cavitation. 
The cavitation is departure from evaporation. The evaporation, 
in definition, is performed by temperature changing but 
cavitation is performed by pressure changing. Cavitation can 
be observed in a wide variety of propulsion and power systems 
like propellers, pumps, nozzles and injectors [1, 2]. 

Cavitation is categorized by a dimensionless number called 
cavitation number, where it depends on the saturated pressure, 
the flow reference pressure, density and velocity, respectively. 
The cavitation number is defined as follows: 


P 7 Post 


o = ———— 
1. (1) 
7 PaUn 


Usually cavitation formation of a flow is categorized based 
on the cavitation number of the flow. Therefore, experimental 
observations are classified based on the cavitation number of the 
flow [3]. Undesirable aspects of cavitation are erosion, structural 
damages, noise and power loss in addition to beneficial features 
such as drag reduction and the effects of cavitation in water 
jet washing systems. The drag reduction observed on bodies 
surrounded fully or partially with cavity strongly encourages 
one to research on properties of cavitating flows. 

In spite of its longstanding practical importance and rich 
physics, cavitating flow continues to be a topic of significant 
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challenge to the computational community. The simultaneous 
presence of interfacial dynamics, multiple timescales, and phase 
change complicates the fluid physics and requires substantial 
modeling efforts. Numerous modeling strategies have been 
proposed in the literature, ranging from Rayleigh-Plesset type 
of bubble formulation, Kubota et al. [4], which separates the 
liquid and vapor region based on the force balance notion, 
to homogeneous fluid approach, Senocak and Shyy [5], 
which treats the cavity as a region consisting of continuous 
composition of liquid and vapor phases. 

In this study the bubble dynamic model which is based on 
the Rayleigh equation is used for simulating phase change. 
A finite-volume approach written in body fitted curvilinear 
coordinates, on collocated grids in 2D and 3D domains is 
used for the numerical descretization. For the present study, 
the transport equation-based model, TEM, is implemented into 
the solver and related modifications, regarding the convection 
schemes and the PISO algorithm, have been made for time- 
dependent computations. 

The main numerical results of this study are from simulation 
of sheet cavitation around NACA66(MOD) hydrofoil and 
supercavitation around a flat plate. Unsteady simulation of 
cavitating flows around a 2d flat plate oriented normal to flow 
is performed from range of subcavitation to supercavitation. 
Apart from their theorical interest (in stability analysis for 
instant), unsteady supercavitating flows were mainly considered 
in the past for their applications to hydrofoils under transient or 
periodic conditions. On the whole, the basis of the modeling of 
unsteady 2d supercavities is the same as for steady flows. 


The main analysis methods are: 

1. analytical, non-linear methods: Von Karman [6], Woods 
[7], Wu [8] 

2. analytical, linearized methods: Timmam [9], Geurst [10], 
Tuln et al. [11] 

3. numerical, non-linear methods, which are currently 
employed and usually use the scheme initially proposed 
by Plesset et al. [12]. A survey of numerical techniques for 
unsteady cavity flow modeling was given by Kinnas [13]. 


Due to the difficulty of conducting tests under unsteady 
conditions, experimental data in this field are rather rare. 
Therefore, the published experimental results usually consist 
of time-averaged properties [14]. 

The simulated flows around the flat plate consist of 
unsteady cavitating flows from two-phase vortex shedding 
to fully supercavitating flows. Based on these simulations, 
the pressure distribution, cavity region, cavity characters 
as its length and width are presented. Beside these, the 
interaction between vapor phase and flow around the flat 
plate, formation of supercavitation and its effects are analyzed. 
Finally, the obtained results are compared with available 
experimental results to demonstrate the accuracy of the current 
simulation. 


GOVERNING EQUATIONS 


The equations governing the flow ofa compressible fluid are 
the continuity equation, the momentum equations, the energy 
equation, and finally the state equation. This set of non-linear, 
coupled equations is solved for the unknown’s parameters p, 
U, T and P. In index notation form, these equations may be 
written as: 


OP m + OP mi; =0 (2) 
ôt Ox; 
OP ni ) OP nuiu j) S 
ot Ox, 
j 6) 
ôP) of (õu â; 
r ara Gl) oe ma | | 
A, A A, 
CT) 4 O nu jT) — 
a» a! 
ot OX; (4) 
1 | ð oT al oP aa as 
| KK BT Ht | OG 
Cp | Ox Ox; Ot OX; 


Where ¢ denotes the dissipation term in energy equation 
and f is the thermal expansion coefficient which is equal to 1/T 
for an ideal gas. In addition to the above differential equations, 
an auxiliary equation of state relating density to pressure and 
temperature [p = f(P,T)] is needed. In many practical problems 
related to cavitation phenomena, the change in temperatures 
is negligible. Therefore the simulation of cavitation in 
isothermal condition has not any effect on final results, and 
it is unnecessary to solve the energy equation. Therefore the 
pressure-density coupling is complex and requires special 
attention [15]. 

In this study, instead of a state equation, the TEM of 
vapor is employed. With this equation, pressure and density 


are connected implicitly via a phase change source term. To 
simulate phase change between vapor and liquid, a term S, , is 
added to the right side of vapor volume fraction equation. The 
vapor volume fraction equation with the phase change source 
term in its right side is presented in equation (5). 


Oty oes )_ S, (5) 
ot Ox; 


The source term of the vapor volume fraction equation 
presents the rate of phase change between vapor and liquid 
phases. The utilized phase change model is presented in the 
part of cavitation model. 


The Homogenous model 


With determination of the volume of fraction, the local 
properties of fluid can be achieved based on the single state of 
each phase. This method is named The Homogenous model. In 
two phase flows, the mixture density and mixture viscosity are 
defined as follows based on the vapor volume fraction: 


Pm = AyPy + (l Ay )Pi (6) 
Um =O Hy +d- ay Ye 


Cavitation model 


In TEM approach, numerical models of cavitation differ 
in cavitation source term, S, . The cavitation source term 
defines vapor net mass generation that contains effects of 
vapor production and destruction. In this study we consider 
the bubble dynamics method as the phase change model [16, 
17]. Therefore S, may be written as: 


. n 
Sa = sign(P... = P) = : 
1+M95 nR 
The average nucleus per liquid volume is considered as 
n, = 10°. Other properties such as the minimum radius of bubble 
can be calculated based on the value of the n, [18]. 


DISCRETIZATION METHOD 


A finite volume method is used to discretize the governing 
equations. The details related to the finite volume discretization 
methods of the Navier Stokes equations and the conservation 
of mass equation can be found in different references [19]. 


Equations descretization 


In the descretized form of Navier-Stokes equations, there 
are three major terms: the unsteadiness, convection and 
diffusion terms. The discretization of the diffusion flux does 
not require any special consideration and the method adopted 
here is the second order estimation. The discretization of the 
convection flux is, however, problematic and requires special 
attention. In this study, the First-Order Upwind (FOU) scheme 
is used to calculate the faces values. However, for increasing 
the solution stability, high order methods could be utilized such 
as Jassak method [20]. For the representation of the unsteady 
term, the time derivative is approximated using the First-Order 
Euler-implicit formulation [16, 20]. 

To discretize the volume of fraction transport equation, it 
is necessary to compute the values on the computational cells 
interfaces accurately. In high speed flows and for capturing 
shock, it is necessary to use high order methods such as 
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HRIC [21, 22]. In low speed flows, the first order Upwind 
method can be used to calculate the volume fraction values on 
computational cells faces [15]. With the homogeneous approach 
and volume fraction values on faces, the density on the cell 
interface can be calculated. 

The velocity fluxes on computational cells faces are 
calculated by using momentum interpolation. The momentum 
interpolation avoids pressure oscillations in the solution 
procedure [21]. 


Overall Solution Procedure 


Selecting a suitable time step for unsteady simulation of 
cavitating flows is very important. The selected time step should 
be proper for convection in the vapor transport equation and 
the continuity equation, as well as Navier-Stokes equations. 
If the time step is selected inappropriately, final results may 
be wrong or the solution procedure may diverge. In some 
researches, the time step is considered in accordance with the 
non-dimensional time of the flow. The non-dimensional time 
is obtained by dividing the length scale to reference velocity of 
the flow [23]. In another approach, the time step is calculated 
in the beginning of each time step by considering the courant 
number, CFL [18]. In this study, the CFL parameter is used to 
calculate the time step. 

After calculating suitable time step, it is possible to start 
the solution procedure. For each time step, first, the vapor 
fraction transport equation is solved and a new vapor fraction 
distribution is obtained. Consequently the values of the 
mixture density and viscosity are updated. Based on these new 
values, the Navier-Stoks equations and the pressure correction 
equations are solved until a convergence criterion is reached. 
Then, the whole procedure is repeated within the next time 
step. In this study, for solving velocity-pressure coupling, the 
non-conservative PISO algorithm is used. 

For the numerical simulation of cavitating flows, the 
pressure level usually is defined by a pressure boundary 
condition at the outlet of the computational domain, and 
velocity set as an inlet velocity boundary. 


NUMERICAL RESULTS 


The numerical model was implemented in the CFD-Code 
developed at the Marine Engineering Lab. of Sharif University 
of Technology. The accuracy of this software are evaluated 
through different numerical simulations. 

The flow is assumed isothermal and fluid properties are 
supposed to be constant at a given temperature for the entire 
flow domain. For all simulations presented, cold water at 
a constant temperature T = 293.2 K with 10° nuclei per m? 
water having a minimal nuclei radius of 30 microns is assumed 
to match the experimental conditions. The saturated pressure, 
P „ is set to a constant value of 2340 Pa. 


NACA66(MOD) Hydrofoil 


Effects of leading edge cavitation on NACA66(MOD) were 
experimentally investigated by Shen and Dimotakis [24]. A2D 
NACA66(MOD) airfoil with camber ratio of 0.02, mean line 
of 0.8 and thickness ratio of 0.09 is used in this simulation. 
The implied boundary conditions and non-orthogonal meshes 
are shown in Fig. 1. 

The published experimental results contain the distribution 
of static pressure on hydrofoil surface at different angles of 
attack and Reynolds numbers. In this study, Leading Edge 
Cavitation simulations were performed at Re = 2 x 10°, an 
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angle of attack of 4 deg, inlet velocity 2.5 m/s and cavitation 
number of 0.91; under these conditions, the cavitation is 
confined to the front of the hydrofoil. Calculated Cp values, 
equation (8), on hydrofoil top surface are shown in Fig. 2 
together with experimental data, and good correlation is clear 
between these results. The vapor volume fraction distribution 
is shown in Fig. 3, which shows the cavitation zone on the 


hydrofoil surface. P -P 
C, =7=— 
2 Po U æ% (8) 
a) Wall 
> : 4 =I Pressur 
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Fig 1. Implied boundary conditions (a) and non-orthogonal meshes (b) 
for simulation of cavitation around NACA66(MOD) 
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Fig. 2. Comparison between obtained numerical pressure coefficient 
distribution and Experimental Results[5] around NACA66(MOD), 
Inlet velocity = 2.5 m/s. Outlet cavitation number = 0.91 


Fig. 3. Vapor volume fraction distribution around NACA66(MOD), 
Inlet velocity = 2.5 m/s. Outlet cavitation number = 0.91 
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Fig. 4. Implied computational domain, boundary conditions and coordinate system for simulation of cavitation around the flat plate 
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Fig. 5. Vapor volume fraction and pressure distribution around the flat plate Inlet velocity = 5 m/s. Outlet cavitation number = 1.0 
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Flat Plate 


Geometry and Boundary Conditions 


The used computational domain and geometry are presented 
in Fig. 4. The plate has height H and thickness 0.1H where H 
is equal to 0.1 m. The boundaries of upstream and downstream 
are located at 10H and 30H. The upper and lower boundaries 
are located at 5H from the center of plate. The constant velocity 
inlet and constant pressure outlet are considered as inlet and 
outlet boundaries. The main results are obtained by considering 
outlet cavitation number equal to 1.0 and 1.25, and velocity 
inlet equal to 5m/s. 


Unsteady behavior of cavitation 


For considering unsteady behavior of fully cavitating 
flows behind the flat plate, the formation of cavitation and its 
oscillations are presented in Fig. 5. The considered cavitation 
number is equal to 1.0 at the outlet boundary. As presented 
in Fig. 5, supercavitation is formed behind the flat plate in 
these conditions. At the end of this cavity, the two phase 
vortex shedding occurs, and vapor is separated from cavity 
by vortexes and moved to the down stream. The cycle of this 
separation occurs in 17.5 ms. Therefore, the frequency is equal 
to 57.14 Htz. Along with the vapor contours, pressure contours 
are presented in Fig. 5, which shows the same oscillating 
behavior. 


Cavity Dimensions 


One of the most important characters of supercavitating 
flows is their cavity dimensions. The cavity dimensions are 
usually normalized by the height of the flat plate. The vapor 
iso-surfaces are presented in Fig. 6. In this figure, one snapshot 
of periodic behaviors of each cavitation number is presented. 
By using these vapor iso-surfaces, average cavity length and 
width can be estimated. Besides, non-dimensional cavity 
length and width are presented in Fig. 7 and Fig. 8. At the 


small values of cavitation number, most of the flow domain 
contains vapor phase. In these situations, the re-entrant jet and 
vapor volume fraction distributions inside the cavity are much 
more complicated than sub-cavitation flows. The obtained 
results for cavity dimensions are compared with experimental 
results reported by Waid [25]. One of the basic objectives of 
this paper is to improve simulation of supercavitating flows. 
For this reason, the obtained results are also compared with 
other numerical simulations [26]. By this comparison, it is clear 
that the implied method improves accuracy of supercavitation 
simulations considerably. 
6 
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Fig. 7. Cavity width defined by a = 0.5 for several cavitation number 
Experimental data by Waid [25], Shin et al. numerical results[26] 


Drag Forces 


Drag force variation is another important character of 
cavitating flows. In the flow around of the flat plate, drag force is 
mainly pressure drag force. The average pressure distributions 
around the flat plate is presented in Fig. 9. By integrating these 
pressure distributions, average drag forces can be obtained. 
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Fig. 6. Vapor volume fraction Iso-surface and cavity dimensions around the flat plate Inlet velocity = 5 m/s, Outlet cavitation number = 1.0, 1.25 
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A Experiment by Waid 
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Fig. 8. Cavity length defined by a = 0.5 for several cavitation number 
Experimental data by Waid [25], Shin et al. Numerical results[26] 
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Fig. 9. Average pressure coefficient distribution of the 2d flat plate for 
several cavitation number 


In the flow around the flat plate, the drag force consists of 
the front force implied in the front of plate, and the back force. 
In the Fig. 10, the unsteady behavior of these two drag forces, 
and total drag force are presented for one cycle. From this 
figure, it is obvious that variations of these forces are similar. 
The average forces for the front and back of the plate, and 
also the total drag force can be obtained from this figure. By 
using the average of the total drag force, part (c) of Fig. 10, 
and the reference velocity and density values, the average 
drag coefficient can be calculated. In the Fig. 11, variations 
of average drag coefficient against cavitation number are 
presented. In this figure, the free stream prediction [12], Shin 
et al. numerical results [26] and experimental results by Waid 
[25] are presented together with present numerical results. 

Beside the undesirable aspects of cavitating flows such 
as noise and erosion, the developed cavity can reduce drag 
forces, especially in supercavitating flows. This beneficial 
aspect of cavitating flows can help to increase the velocity 
of vehicles in constant power. When speed of vehicle isn’t 
enough to create supercavitation naturally, Ventilation can be 
used to create or to enhance a supercavity called ventilated 
cavitation. In the Fig. 11, it is presented that by development 
of cavitation and reduction of cavitation numbers of the flow, 
drag force reduces. The summery of obtained numerical 
results for simulation of cavitation around the flat plate are 
presented in Tab. 1. 
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Fig. 10. Unsteady behavior of the front (a), back (b) and total drag force (c) 
of the 2d flat plate Inlet velocity = 5 m/s, Outlet cavitation number = 1.0 
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Fig. 11. Average drag coefficient of the 2d flat plate for several cavitation 
number Experimental data by Waid [25], Shin et al. Numerical results[26], 
Free Stream Theory[12] 
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Tab. 1. The summery of obtained results for simulation of cavitation around the flat plate 


| oo | Re | um | oH | Ave. Drag Force (N) | Ave. Drag Coefficient 
1.0 6.2 2.2 1:75 


2193 


2888 


Tab. 2. The comparison between cavitating and non-cavitating simulation results 


Cavitation Simulation Results 


Non-Cavitation Simulation Results 


Ave. Drag Force (N) 


It seems necessary to estimate gain accuracy via simulation 


of cavitation in a flow. For this reason, the simulated cavitating 
flows are compared with results of non-cavitating simulations. 
This comparison shows simulation of cavitation can provide 
more accurate results for drag force, periodic behavior, and 
pressure distributions. The comparison between caviatiting 
numerical results and non-cavitating numerical results are 
presented in Tab. 2. 


CONCLUSIONS 


Unsteady simulations have been performed for 2D 
cavitating flows. For phase change modeling, the Bubble 
dynamics cavitation model is utilized, which is presented 
as the source terms in the volume of fraction equation and 
the continuity equation. The non-conservative PISO method 
is used to solve coupling between the continuity and N.S. 
equations. Unsteady behavior of cavitation around a flat 
plate oriented normal to flow direction, and sheet cavitation 
around a NACA66(MOD) hydrofoil are performed to 
demonstrate the accuracy of the selected algorithm. 


Based on the implied phase change model, the supercavitating 
flow behind the flat plate has been simulated in a good 
agreement with the experimental results. The size of 
the vapor region and the separation of the cavitation are 
clearly predicted by the developed CFD-code. These 
predictions provide information that will be helpful for 
understanding behaviors of unsteady cavity flows such as 
cavitation occurrence and development. The effects of vapor 
generation on drag reduction are investigated and presented. 
It is presented that by growth of cavity and reduction of 
cavitation number, the drag force reduces. 


The results of the unsteady simulation show the development 
of a re-entrant jet and two-phase vortex shedding. The 
predicted cavity dimensions, drag coefficients and pressure 
distributions show a very good agreement with the 
experiments. 


NOMENCLATURE 

C  — hydrofoil chord 
CFL — courant number 

C, — pressure coefficient 
g;  — gravity 

H — plate height 

K 


— termal expansion coefficient 
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Ave. Drag Force (N) | Ave. Drag Coefficient 


— cavity average length 

— average nucleus per liquid volume 
— pressure 

— saturation pressure 

— energy source term 

— gas constant (eq. ()) 

— nucleus radius 

— reynolds number 

— vapor phase change term 

— time, mean flow time scale 

— temperature 

— cavity maximum thickness 

— cartesian velocity components 
— velocity 

— cartesian coordinates 


— volume fraction 

— termal expansion coefficient 
— cavitation parameter 

— energy dissipation term 

— viscosity 

— density 


subscripts, superscripts 


œ — free stream 

m — mixture 

l — liquid phase 

v _— vapor phase 

sat — saturation condition 
i,j — coordinate indices 
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ABSTRACT 


The article presents a concept of a combined large power ship propulsion system consisting of the leading 
Diesel main engine, associated with a power gas turbine and the steam turbine system which utilise the 
energy contained in the exhaust gas leaving the Diesel engine. In the examined variant of the combined 
system the power turbine is fed in series with the exhaust gas. A calculation algorithm is given along with 
the results of calculations of particular subsystems of: the turbocharging system, the power gas turbine, and 
the steam turbine cycle. Assumptions and limits adopted in the calculations are presented. Selected system 
parameters were confronted with the results of experimental investigations available in the literature. The 
power optimisation of the entire combined ship power plant was only performed taking into account the 
thermodynamic point of view, leaving aside technical and economic aspects. Numerical calculations were 
performed for the 52 MW low-speed Diesel engine produced by Wärtsilä. 


Keywords: Ship power plants; combined systems; Diesel engine; gas turbine; steam turbine 


INTRODUCTION 


The low-speed main engine used for driving a ship has the 
highest efficiency of all heat engines. What is more, it burns 
the cheapest residual fuel. However, there are some problems 
with keeping the emission of harmful substances by this engine 
to the atmosphere at a permissible level. 

The article analyses a variant in which the turbocharger 
and the power turbine are connected in series. An option was 
also taken into account in which the exhaust gas was utilised 
in the waste heat boiler working for the steam turbine. 

The analysed combined power plant consists of a Diesel 
engine 9RTA-96C, identical to that discussed in [2, 3], the 
power turbine, fed in series, and the steam system with 
a steam turbine. Additional power of the propulsion system, 
obtained both from the power turbine and the steam turbine, 
is transmitted to the ship propeller shaft via a mechanical gear. 
Behind the turbocharger, a power gas turbine is installed in 
series, from which the exhaust gas flows through the waste heat 
boiler producing the steam to drive the steam turbine. Steam 
consumption for ship’s own needs was taken into account in 
the calculations. 

For the assumed efficiency of the turbocharger and the 
power turbine, the power of the gas turbine was calculated, 
while the steam turbine cycle was optimised with respect to 
the maximum power of the adopted two-pressure waste-heat 
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boiler. The optimisation of the steam turbine system took into 
account limits resulting from practical solutions used in power 
transmission systems. 


CONCEPT OF A COMBINED SYSTEM 
— SERIES FEEDING OF POWER TURBINE 


The combined cycle includes an additional power turbine 
and a waste heat boiler which produces the steam for ship’s 
own needs and for an additional steam turbine. The use of 
the power gas turbine is possible due to high efficiency of the 
turbochargers, which makes it possible to increase the pressure 
of the gas behind the turbocharger turbine. This pressure is 
higher than the pressure of the exhaust gas in front of the waste 
heat boiler, which results from pressure losses connected with 
the flow of the exhaust gases through the boiler. The pressure 
difference between these two pressures provides opportunities 
for the use of an additional gas turbine, bearing the name of 
the power turbine, in the system. 

This solution leads to the increase of power of the entire 
power transmission system, as the volume of fuel delivered to 
the Diesel engine does not change, while both the power turbine 
and the steam turbine deliver additional power. 

It results from operating low-load tests of conventional 
power transmission systems with Diesel engines with 
turbocharging that the amount of the exhaust gas in the 


Fig. 1. Combined ship power transmission system 


turbocharger is not sufficient to secure stable engine operation. 
In that case additional air blowers, situated at engine inlet, are 
started. That is why the proposed combined power transmission 
system should have power turbine bypass pipes (Fig. 1) which 
direct the exhaust gas from the turbocharger to the waste heat 
boiler during low-load intervals. This system will be controlled 
automatically by a valve for which the control signal will be 
the pressure of the turbocharging air, the engine torque, or 
the angular speed of the propeller shaft. This control system 
also seems to provide opportunities for the operation of the 
ship power transmission system without the power turbine. 
A possibility to switch off the power turbine results in the 
increase of the turbocharger‘s power output, which further 
leads to the increase of both the turbocharging pressure and 
the volume of the turbocharging air, all this increasing the 
ability to manoeuvre the engine, in particular fast reaching the 
maximum power at important moments of power transmission 
system operation points. Also in cases of large power of the 
main (piston) engine, increasing the power of the power turbine 
may lead to the decrease of the power of the turbocharger, thus 
increasing the power of the entire power transmission system. 


OPERATING PARAMETERS OF THE MAIN 
ENGINE TURBOCHARGER 


The engine 9RTA-96C has constant-pressure turbocharging 
[5]. Three PCL-type turbochargers were produced by ABB. 
From the exhaust gas collector the exhaust gas, having the 
temperature t „p and the pressure P p p flows directly to the 
turbochargers. In the catalogue the producer [5] does not present 
all parameters concerning turbocharging. In further calculations 
the temperature t , p and the pressure P p p of the exhaust gas 
in the exhaust collector were assumed based on the analysis 
reported in [2]. 

It was assumed in turbocharging calculations for the 
combined system that the power needed for driving turbocharger 
turbines is sufficiently large if the pressure behind the 
turbocharger turbine is higher than the inlet pressure of the 
exhaust gas in the waste heat boiler. 


The pressure of the exhaust gas behind the turbocharger is 
calculated from its power balance: 


l 
Nr = 
T ] -l 
m c t LS a 
= z air Ya, bar (E) 
Nre Mre Cg tinet TC 
Pexh_ D 
Pexh Tc 5—77 (1) 
Tr 
where: 
To= Pard Pis is the compression ratio in the turbocharger 
compressor 
T= Pon r Ponp 1S the expansion ratio in the turbocharger 
turbine. 


The efficiency of the turbocharger was defined as: 
Nrc = Nc ` Nr ` Tii (2) 


The compressor efficiencies n, and n, in the above formula 
were determined in the following way: the compressor 
efficiency was calculated from the line of cooperation of 
the engine RTA-96C with the ABB compressor. Changes of 
turbocharger turbine efficiency were calculated from MAN 
B&W experimental tests [4]. It was assumed that the efficiency 
of the ABB turbine changes in the similar way. This relation 
made the basis for determining the efficiency of the ABB 
turbocharger turbine cooperating with the engine 9RTA-96C 
as a function of its load. The results are shown in Fig. 2. 

The temperature of the exhaust gas leaving the turbocharger 
turbine was calculated from the transformed formula which 
defines the turbine efficiency as: 


ai —273.15 
Kg—l 


Tre “e (3) 


texh to= (texn p+ 273.15): l-nr|l- 
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Fig. 2. Characteristics of turbocharger parameters as function of main 
engine load: £ - turbocharger outlet temperature in a standard engine 
(producer 5 data), x* - expansion rate in the turbocharger turbine 
(producer 5s data) 


It results from the performed calculations that when the load 
is larger than 65%, the operation of the power turbine in the 
system is possible, Fig. 2. For the presently used turbochargers 
the turbocharger turbine power is reached at low compression. 
The excess pressure of the exhaust gas can be expanded in the 
additional power turbine. 


POWER TURBINE CALCULATIONS 


In the adopted concept of the combined power transmission 
system, Fig. 1, the power turbine is situated behind the 
turbocharger. The exhaust gas from the turbocharger flows to 
the power turbine. The entire mass flow rate of the exhaust gas 
leaving the main engine expands in the power turbine, thus 
delivering extra power to the ship power transmission system. 

The pressure behind the power turbine depends on the losses 
taking place during the flow through the waste heat boiler and 
outlet silencers. Following common practice, it was assumed in 
further calculations that the pressure behind the power turbine 
is by 3% higher than the barometric pressure, i.e. 


Pexh PT E 1.03 Poar (4) 


The expansion ratio in the power turbine can be 


calculated as: p 
inlet PT 
ter 5 (5) 
Pexh PT 


The power of the power turbine was calculated from the 
following formula 


Npr = Mer’ Mpr' Cg‘ tinter pr |1— 


Ko 
(Ter) * 
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It was assumed in the calculations that the power turbine 
efficiency npp depends on the expansion ratio, according to the 
data presented by MAN B&W for the turbocharger turbine [1]. 
This assumption can be considered correct, as the producers 
make use of gas turbines installed in turbochargers as power 
turbines. 

The temperature of the exhaust gas leaving the power 
turbine was calculated from the following relation: 


—1_}-27315 
Kg—l 


Tpr “eg (7) 


texh pc=(tintet prt273.1 5): at (| 


The power turbine parameters were calculated using 
a computer programme (the source code is included as Appendix) 
taking into account the above assumptions. Fig. 3 presents the 
calculated values of the selected power turbine parameters as 
a function of the internal combustion engine load. 

The use of the power turbine in the system resulted in 
power increase of the entire system from 3 to 9 %, compared 
to the standard engine power. The power of the power turbine 
increases with the increase of power of the internal combustion 
engine. It results from the performed calculations that when 
the load is smaller than 65 % of the internal combustion engine 
power, the power turbine cannot be used as its power nears 
zero. In this case the power turbine is to be disconnected from 
the system. This conclusion is in accordance with the results 
of turbocharging process calculations. 


Pintet_PT 


p [bar] 


55 65 75 85 95 105 115 
Np [%] 


Fig. 3. Selected power turbine parameters vs. main engine load 


STEAM CYCLE OPTIMISATION 
CALCULATIONS 


In the analysed combined power transmission system, 
Fig. 1, the exhaust gas leaving the power turbine is directed 
to the waste heat boiler, where it produces steam used to meet 
ship’s own needs and supply the steam turbine. The steam 
turbine transmits the power, via a gear, to the propeller shaft. 
The system of the steam turbine cycle was adopted following 
the description given in [2]. In the steam system a two-pressure 
waste heat boiler was used. The calculations were preformed 
using the computer programme [2] to search the maximum 
power of the steam turbine. The steam cycle was optimised 
for two main engine loads, i.e. for the power corresponding to 
100% and 90% of the thermal load of the piston engine. 

Table 1 collects the results of calculations for two internal 
combustion engine load regimes which corresponded to the 
maximum available power of the steam turbine. It results from 
the analysis of the steam cycle calculations performed for 
the examined combined power transmission system that the 
maximum available power of the steam turbine equals 7.25 


or 7.96 % of the internal combustion engine load, depending 
on the load. 

In practical application, the steam turbine system needs 
some limits to be placed on its parameters as a result of the 
adopted designing, constructional and operating assumptions. 
For the adopted limits [2] the steam power transmission system 
was optimised for the second time to find a solution for which 
the power of the steam turbine reaches its maximum. Tab. 1 
presents the parameters of this system, obtained taking into 
account the adopted limits. Due to the limits in the steam 
turbine cycle we cannot utilise the entire power resulting form 
the calculations performed for the cycle without limits. The 
presented analysis of the steam cycle reveals that the additional 
use of the steam turbine in the combined system increases the 
power output of the system by about 6.5 + 7.8 % and decreases 
the fuel consumption. 


THERMODYNAMIC ANALYSIS OF THE 
COMBINED POWER TRANSMISSION 
SYSTEM 


The application of the combined system in a classical 
ship power transmission system with a low-speed Diesel 
engine leads to the increase of power and to the decrease of 
the specific fuel consumption. Tab. 2 collects the results of 
calculations of particular components of the combined system. 
It results from the performed analysis that the sole use of 
a power turbine increases the power output of the system by 
5.5 + 7 %, depending on the Diesel engine load. At the same 
time the steam system alone (without the power turbine) 
increases the power output of the plant by 6.5 + 7.5 %. In the 


complete combined system composed of the piston engine as 
a leading engine, the power turbine, and the steam turbine, the 
power output of the system increases by 12 + 15 %. This power 
output increase is not generated by extra fuel delivered to the 
system but is connected with deep utilisation of the exhaust 
gas leaving the piston engine. This system also makes it 
possible to reduce the specific fuel consumption by 11 + 13 % 
compared to the traditional power plant with an internal 
combustion engine. 


COMPARING VARIANTS OF SHIP POWER 
TRANSMISSION SYSTEMS 


Tabs 3, 4 and 5 show selected parameters of the calculated 
combined systems, which were taken from the above presented 
variants of the combined ship power transmission system 
composed of the main engine, and the power and steam turbine 
systems. 

In the steam turbine cycle of the combined system, the 
cycle parameters differ little between each other in the 
analysed variants. The live steam temperature differs by 
3°C in one case, while the live steam pressure also does not 
differ much. The pressure in the first stage of the boiler and 
in the degasifier, as well as the temperature of the exhaust 
gas leaving the power transmission system are the same for 
all variants, Tab. 3. The largest power output of the steam 
turbine cycle is reached in the combined system variant with 
the steam turbine and power turbine fed in series with the 
exhaust gas, Tab. 4. Differences in steam turbine power output 
for particular variants for 90% of the main Diesel engine load 
are approximately equal to 4%. 


Tab. 1. Optimum steam cycle parameters in a combined cycle consisting of Diesel engine - power turbine (series feeding) - steam turbine 


without limits 


154 


82 


153 


82 


limits 


121 


Tab. 2. Power output, efficiency and specific fuel consumption of the combined ship power transmission system 
with the two-stroke engine 9RTA 96C, power turbine (fed in series) and steam turbine 


ecombi 
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Tab. 3. Comparing steam cycle parameters for load equal to 90% of main engine power 


Combined system 


Steam turbine 


Power turbine - fed in parallel 


Power turbine - fed in series 


Tab. 4. Comparing power outputs for different variants of the combined system — with power output of the standard main engine as the reference 


Main engine — 


Main engine — power turbine — steam turbine 


Main engine power output steam turbine [2] 


Steam 


Feeding in parallel [3] Feeding in series 


turbine 


Power 


turbine 


For combined system variants with the power turbine, the 
largest power is obtained in the variant with the power turbine 
fed in series, Tab. 4. In this cases the power output is higher by 
about 16.3 % compared to the power turbine from the variant 
of parallel turbine feeding with the exhaust gas for 90% of 
main engine load. 

In the analysed variants the power output of the steam 
turbine is larger than that of the power turbine by 6.1 - 28.8 %, 
depending on the system variant and main engine load, 
Tab. 4. 

Depending on the variant and main engine load, the 
application of the combined ship power transmission system 
makes it possible to increase the power output of the power 
plant by 6.9 — 14.6 % as compared to the conventional power 


plant, without delivering extra fuel, Tab. 5. The additional 
power output of the system is obtained by utilisation of the 
energy in the exhaust gas leaving the piston engine. This way 
the combined system reduces the specific fuel consumption by 
6.4 — 12.8 % compared to the conventional power plant. 

Depending on the adopted solution, the application of the 
combined power plant makes it possible to reach the assumed 
power output of the power transmission system at smaller main 
engine loads, thus reducing the 24-hour fuel consumption of 
the ship power plant. Tab. 6 presents, for the assumed power 
output of the ship power plant, the main engine power output 
and the reduction of the 24-hour fuel consumption, in relation 
to the conventional power plant, as a function of the adopted 
variant of the combined power plant. 


Tab. 5. Comparing different combined system variants 


Main 
engine load 


Main 
engine 


System power 


Combined power transmission system 


Steam Power turbine fed 


turbine [2] 


In parallel [3] In series 


output Lewd 


AN 


combi D 


[%] 


be. /be [g/kWh] 


combi D 
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Efficiency q [%] 
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Tab. 6. Reduction of 24-hour fuel consumption of the 49500 kW ship power plant for particular combined system variants 


Conventional power plant 


Combined system 


D & PT & ST - feeding: 
In parallel [3] 


D & ST [2 
2] In series 


Main engine load 


24-hour fuel consumption 


Fuel consumption reduction 
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FINAL CONCLUSIONS FROM THE 
ANALYSIS OF COMBINED SHIP POWER 
TRANSMISSION SYSTEM VARIANTS 


It results from the analysis of the optimised combined 

ship power transmission system that the highest efficiency 

is obtained when particular power transmission system 

components are used: 

— diesel engine with maximum efficiency 

— turbocharger. Most highly efficient turbocharger is to be 
used, to provide opportunities for reducing the exhaust 
gas pressure at its outlet, and, consequently, increasing 
the power of the power turbine 

— power turbine. High efficiency power turbine is required 
to increase its power output 

— steam turbine cycle. The waste heat boiler is expected 
to reveal low exhaust gas flow losses (they decrease the 
expansion end pressure in the power turbine) and small 
temperature concentration in boiler evaporators (pitch 
point). Observed is the high impact of sulphur content 
in the fuel on the permissible temperature of the exhaust 
gas and the low temperature limit for the supply water. 
The least possible number of heat exchangers (none, as 
an ideal case) is to be used in the steam turbine cycle. 
The optimum parameters of this cycle also depend on 
the piston engine load. 


The power transmission system of the combined ship power 

plant of a big transoceanic container ship, consisting of 

the main Diesel engine, the power turbine and the steam 
turbine, makes it possible to: 

— increase the power plant power output by 6.9 + 14.6 % 
compared to the conventional power plant, without 
delivering extra fuel 

— increase the power efficiency of the power plant from 
53 % to about 56 % (which, so far, was impossible 
to reach by the internal combustion engine alone) 
depending on the applied variant of the combined 
power plant. Reducing the content of sulphur in the fuel 
provides opportunities for reaching higher efficiencies 
resulting from larger power of the steam turbine — lower 
temperature of the exhaust gas leaving the waste heat 
boiler 

— reduce the specific fuel consumption by 6.4 + 13 %, 
without delivering extra fuel, as a consequence of 
the increased power output compared to that of the 
conventional power plant 

— use the internal combustion engine having lower power 
output, due to the extra power obtained in the combined 
system, or to reach 100% of power output of the power 
plant when the Diesel engine load is at an approximate 
level of 90%. 


In the article, only thermodynamic analysis of the power 
plant is presented, without additional technical and 
economic analyses, which could fully justify the application 
of this power plant system for driving a large container 
ship. 


NOMENCLATURE 

b, - specific fuel oil consumption 

CC, - specific heat of exhaust gas and air 

i - specific enthalpy 

m - mass flow rate 

N - power 

p - pressure 

T,t - temperature 

Wu - caloric value of fuel oil 

x - steam dryness ratio 

n - efficiency 

KK, - isentropic exponent of exhaust gas and air 

Indices 

bar - barometric conditions 

B - boiler 

combi - combined system 

D - marine diesel engine, supercharging 

f - fuel 

inlet - inlet duct 

k - parameters in a condenser 

o - live steam, calculation point 

air - air 

ss - ship living purposes 

T - compression ratio in a compressor, expansion ratio in 
a turbine 

C - compressor 

g - exhaust gas 

T - turbine 

TC - turbocharger 

PT - power turbine 

ST - steam turbine 

exh - outlet duct 

FW - water supplying a waste-heat boiler 
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Assessing diagnostic applicability of heat release 
characteristics determined based on ship engine 
indicator diagrams 


Stanistaw Polanowski, Prof. 
Polish Naval University 


ABSTRACT 


Inorder to determine indicator diagrams-based heat release characteristics, a single- zone 
model of net heat release was used for perfect gas. It was proved that when a constant 


N 


value for isentropic exponent is assumed, the error in determining the characteristics 
can be limited to 1% at the nominal load. The effect of errors in determining the position 
of piston TDC, as well as that of gas passages and indicator valves on the calculated 
characteristics was evaluated. It was shown that for low-speed engines the effect of gas 
passages in negligible, while for medium-speed engines the characteristics reveal some 


deformations (waves), which are repeatable for an individual cylinder in the examined engine construction. 
The results of the performed investigations and analyses suggest possibility and advisability of the use of 
heat release characteristics in diagnosing ship engines, in particular low-speed machines. 


Keywords: ship piston engines; indicator diagrams; heat release characteristics; diagnostic application 


INTRODUCTION 


Since the production of the combustion pressure analyser 
NK-5 by Autronica, and the analyser Cyldet 1800 by ABB in the 
seventies of the last century, the range of use of the information 
contained in the indicator diagrams has not changed. 

For the time being, the main functions executed by pressure 
analysers include: calculating average indicated pressure, 
determining maximum combustion pressure, recording TDC 
pressure, and the presentation of the determined values in the 
tabularised form and histograms, as well as the presentation and 
comparison of indicator diagrams on a monitor screen. Some 
analysers have also included an option of determining the rate 
of combustion pressure building up. 

Diagnostic inference based on tables of parameters, 
histograms and comparison of indicator diagrams is rather 
unreliable and highly inaccurate. Hence, some monitoring and 
diagnostic systems have included an additional option of pressure 
measurement and analysis in the fuel injection system. 

Perception of the diagnostic information contained in the 
indicator diagram can be extended by making use of heat release 
characteristics determined based on the indicator diagrams. 


SELECTING THE HEAT RELEASE MODEL 


It is Paul Henry Schweitzer [5] who is believed to be the 
author of the first heat release model for the Diesel engine [1, 
4]. Remarkable development of advanced calculation models, 
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observed since the late fifties of the last century, is connected 
with the appearance of possibilities for computer aided data 
processing and computer simulations. 

At present, model investigations of combustion processes 
make use of complicated, multidimensional models of heat 
release processes. For instance, most advanced models include 
simulation programmes KIVA worked out in the Los Alamos 
National Laboratory. These model are characterised by a large 
number of immeasurable quantities and parameters, and their 
use requires high computational potential. 

The simplest and most frequently used model which 
calculates heat release in a Diesel engine based on indicator 
diagrams is believed to be the single-zone model developed 
by Krieger and Borman [1, 4]. It is assumed in this model that 
at each time instant the working medium (charge), having the 
form of uniform mixture of air and combustions products, is 
in the state of thermodynamic equilibrium. 

The first law of thermodynamics in the differential form 
for an open system can be written as: 


dQ- dQ,- dQ. =dQ,=dU+dW (1) 


where: 
dQ, — heat released as a result of fuel combustion 
dQ, — heat of cooling 


dQ, — heat lost as a result of transfer of substances through 
system boundaries 

dQ, — net heat release 

dU — internal energy of the charge 

dW  — work done by the system. 


The heat Q, is the result of fuel flow to the cylinder and gas 
scavenge through piston rings and valves. 

The quantities Q „and Q, cannot be determined in operating 
conditions due to the lack of required measurable data. 
Therefore for diagnostic purposes it is advisable to use net heat 
release characteristics Q 

The formula for heat dQ, for perfect gas can be written as: 


dQ, =(x-1)' (V dp +x p dV) (2) 


where: 

k — adiabatic exponent 
V — gas volume 

p — gas pressure. 


After dividing equation (3) by the displacement volume, 
the formula for the net heat release rate as a function of the 
crankshaft rotation angle takes the form: 


N 
- 4. (x - (vee ae + po” (3) 
da da j 
where: 
v — dimensionless gas volume [3], 
a — crankshaft rotation angle. 


Here, the basic unit for the quantity q' is [J/m?°-° OWK]. 
The net heat release, q, in relation to the stroke volume and 
as a function of the crankshaft rotation angle, is expressed by 


the integral: 
q= fada (4) 


where the beginning of integration 1s assumed at BDC piston 
position. 


ASSESSING THE EFFECT OF CHANGES 
OF THE EXPONENT k ON THE 
ACCURACY OF THE DETERMINED 
CHARACTERISTICS 


For both the air and the exhaust gas, the « values almost 
linearly depend on the gas temperature (Fig. 1). 
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Fig. 1. Adiabatic exponent k vs. temperature T. The characteristic was 


made based on the approximation of tabularised data [7]: 
k,— dry air, k, — exhaust gas 


The effect of changes of the adiabatic exponent «x on q' 
can be taken into account using, for instance, the method of 
successive approximations. For the first approximation we 
assume a constant (average) value of «. In the next calculation 
step the fuel consumption and the volume of the produced 
exhaust gas are assessed based on the heat release determined 
in the precious calculation step. Although it looks simple, this 
algorithm would make the calculation process take much longer 
without delivering any additional diagnostic information. 
The level or the possible error made when we stop at the first 
approximation of the characteristics is to be estimated. 


The temperatures at characteristic points of the working 
cycle for ship engines are within the ranges of : 800-1100 K 
for the end of compression, 1700-2000 K for the end of 
combustion, and 900-1200 K for the end of expansion [6], 
which in total covers the temperature range of 800-2000 K. If 
we assume in the calculations that the value of x is the average 
calculated for the middle temperature within the above range, 
then the errors which we make when determining the maximum 
value of q' do not exceed + 3.5 %, and for q they do not exceed 
+2 %, which was checked on a number of types of ship engines. 
The scale of this error can be even more reduced (below 1 %) 
when we take into account the range of temperatures of the 
working medium which is characteristic for the examined type 
of engine and given load. 


DEFORMATIONS OF CHARACTERISTICS 
INTRODUCED BY GAS PASSAGES AND 
INDICATOR VALVES 


In order to recognise the nature of the effects of the gas 
passages and indicator valves on the heat release characteristics, 
a series of tests were performed on the research rig of the 
medium-speed engine Sulzer 6AL20/24. The arrangement and 
dimensions of the gas passages and the distribution of pressure 
sensors are schematically presented in Fig. 2, preserving 


relevant proportions. 
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Fig. 2. Scheme of indicator passages in the Sulzer engine 6AL20/24, with 
the distribution of pressure sensors: IV — indicator valve, C — Optrand 
pressure sensor in the cylinder V — Kistler pressure sensor 
on the indicator valve 


For technical reasons, the measurements were performed 
using pressure sensors of two different types and made by two 
different producers. The measuring paths were calibrated before 
the measurements and checked after the measurements for 
100% and 50% of the nominal load of the engine. The sensors 
were installed on a special connector fixed to the indicator valve 
in such a way that during the calibration and final check their 
front plates were opposite to each other, at a distance of 5 mm. 
Differences in indications in the assumed range of measurement 
did not exceed + 0.5 % of the maximum value of the measured 
pressure for the nominal load. 

Before calculating the heat release characteristics, the 
indicator diagrams p, and p, (Fig. 4) were smoothed using 
the method of multiple approximation and a third-order 
power polynomial [2]. Without this procedure the heat release 
rate curves would not be readable. The smoothing removed 
(smoothed) the high-frequency measuring noise, Dp, and Dp,,,, 
introduced by the sensors, and the disturbances Dp.» Dp,,, 
generated by gas oscillation in the passages (Fig. 3). 

The above smoothing is necessary to obtain sufficient 
smoothness of the heat release characteristics. However, it 
does not remove the deformation of the pressure time-history 
p, on the indicator valve (Fig. 4), caused by gas compression 
and expansion in the passages, heat transfer, and adding up of 
the remaining errors introduced by the sensors. 

The deformation of the p,, curve, as compared to the pressure 
po in the cylinder, was a source of relevant deformations in 
the heat release characteristics q', and q, , compared to the 
characteristics of q', and q, (Fig. 4). 
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Fig. 3. Disturbances removed (smoothed) from the indicator diagrams after 
their smoothing using the method of multiple moving approximation and the 
third-order power polynomial: Dp, Dp o~ for po (Fig. 4) in the cylinder, 
Dp,» Dp, -for p, (Fig. 4) on the indicator valve 
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Fig. 4. The effect of gas passages and indicator valves (Fig. 2) on the time- 
history of pressure p,,, and the heat release characteristics q', and q, for 
nominal load: Po, q'¢» Ic —for cylinder, p,, q',, q, —for indicator valve, 

W — deformation (waves) observed on q', and q, curves 


Along with smaller peak values, characteristic deformations 
(waves) W can be observed on the q', and q, curves. They 
can also be noticed on the characteristics of medium-speed 
engines of other types (Fig. 7). On the other hand, this type of 
deformation does not exist, in practice, on the characteristics 
recorded on low-speed engines (Fig. 6). 


THE EFFECT OF TDC POSITION 
ESTIMATION ERROR ON THE SHAPES OF 
CHARACTERISTICS 


It is noteworthy that in the case of the examined engine 
(Fig. 2) the deformation of the pressure time-history recorded 
on the indicator valve is accompanied by a pressure signal 
delay, amounting to 2.7? OWK, which was removed by 
shifting the diagrams to their TDC points. In the examined 
case for the measurement on the indicator valve at nominal 
load, the TDC position estimation error, which was equal to 
+ 1° OWK, was a source of errors in determining maximum 
q’ and q values which approximately amounted to + 7 % 
(Fig. 5). 

For engines of different types, the effect of TDC position 
estimation error on the maximum values of q’ and q is 
comparable with its effect on the errors in determining the 
average indicated pressure. The applied method of TDC 
position estimation on indicator diagrams secures high 
accuracy, especially in case of low-speed engines [3]. 


34 POLISH MARITIME RESEARCH, No 3/2009 


q/qomax (A), G/Gomax (B) 


160 180 200 220 240 
a [°OWK] 


Fig. 5. The effect of the TDC position estimation error Aa, on the shapes 
of heat release characteristics q’, and q, (Fig. 4): ye MA ona, We the 
maximum values for the reference characteristics (Aa, = 0) 


ASSESSING THE POTENTIAL FOR 
DIAGNOSTIC APPLICATION OF HEAT 
RELEASE CHARACTERISTICS 


For a newly built and properly adjusted engine, strong 
convergence of q’ and q characteristics is observed for particular 
cylinders (Fig. 6). 
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Fig. 6. Comparing heat release characteristics q' and q for all 5 cylinders 


s =. 


in a newly built low-speed engine 5RTAS52 (sea trials): Tnax’ Vmax 
(referential values) — maximum q' and q values recorded in particular 
cylinders and averaged for the entire engine 


Noteworthy is the coincidence of the waves in particular A 
curves (Fig. 6), irrelevant of the location of the measuring point 
(cylinder). This tendency is also characteristic for other types 
of low-speed engines and for medium-speed engines. 

For an engine after large number of operating hours or 
an inefficient one, the q’ and q characteristics differ between 
each other depending on the technical condition of fuel system 
components (Fig. 7). 

These characteristics (Fig. 7) were determined after the 
running repair and the adjustment of the injection system. 
Remarkable differences in the characteristics are caused by 
advanced wear of injection system components due to a large 
number of operating hours without engine renovation. 

Based on the heat release characteristics we can make 
a judgement about the course of the fuel injection and 
combustion process, and evaluate the accuracy of its adjustment, 
efficiency and wear of particular injection systems. Additional 
information can be obtained by analysing trends of selected 
parameters of the characteristics. 

It is advisable to create a catalogue of symptoms (patterns) 
of typical inefficiencies for engines of certain types. This can 
be done via observation of engines in operation or simple 
simulation tests. 
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Fig. 7. Comparing heat release characteristics q' and q determined for one 
block of cylinders in the 12-cylinder two-stroke medium-speed engine 40DM 
(n = 750 rev/min, nominal load) after large number of operating hours: 
Tnax’ Una referential values) — maximum q' and q values averaged for the 
entire engine, Outlet — exhaust space 

The permeability of exhaust ducts, in particular in two- 
stroke engines which burn heavy fuel, can be diagnosed using 
characteristics parts corresponding to the exhaust space (Fig. 7). 
Some inadequacy of the adopted heat release model for this 


space is of no importance here. 
CONCLUSIONS 


e The net heat release characteristics can be determined 
based on ship engine indicator diagrams, with the accuracy 
sufficient for diagnostic purposes. 

¢ For low-speed engines, deformations of the characteristics 
introduced by passages and indicator valves are negligible. 
For medium-speed engines the deformations of pressure 
time-histories are large but repeatable, which provides 
opportunities for their use for diagnostic purposes. 

e Based on the heat release characteristics one can conclude 
about the range and uniformity of load distribution over 
particular cylinders, condition of the control, correctness 
of work and technical condition of the injection system, as 
well as on the level of exhaust gas resistance. 

e Itis advisable to create a catalogue of inefficiency symptoms, 
which will help ship engineers to detect and recognise 
inefficiencies based on the heat release characteristics. 


NOMENCLATURE 
(repeated symbols which need defining) 


TDC  — top dead centre (of a piston) 

BDC  — bottom dead centre (of a piston) 

Pe — pressure (indicator diagram) measured in the cylinder 

Py — pressure (indicator diagram) measured on the indicator 
valve 

q' — net heat release rate related to the stroke volume, 
[J/m3-° OWK] 

q — net heat release related to the stroke volume, 
[J/m3-° OWK] 

d'o Ie q', q determined for the indication in the cylinder 

darde q', q determined for the indication on the indicator 
valve 

OWK — degree of crankshaft rotation angle (crank angle). 
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The diagnosis of onboard generators 
(alternators) 
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Abstract 


In the paper selected problems related to diagnostics of onboard generators and alternators fitted with 
control systems are discussed. Problems refer to commutator generators and synchronous single- and three- 
phase alternators. Special attention is paid to commutation effects. Results of incorrectness and possibility 
to detect them are discussed. There are also discussed effects associated with changes in a character of 
pulsation, which occur during shortings or insulation clearances in rotor or stator wiring. Possibility of 
diagnosis of generator ’ or alternator 5 parts by means of analysis of pulsation component parameters is 
indicated. In the case of alternators a number of diagnostic methods based on observation of changes in 
shape of voltage or frequency modulation, is discussed. This allows to detect many mechanical or electrical 
faults of generators, alternators or their control systems. 


Keywords: technical diagnostics; frequency modulation 


CHARACTERISTICS OF VOLTAGE 
PULSATION OF DIRECT CURRENT (DC) 
GENERATOR 


In a classical educational approach the DC commutator 
generator is presented in the schematic form, as given in 
Fig. 1, and run of its electromotive force - in Fig. 2. The DC 
commutator generator is consisted of: 

— the motionless stator which can be schematically presented 
in the form of pair of permanent magnets (Fig. 1, where: ,,N” 
— north pole, ,,S” — south pole), producing constant magnetic 
field of B intensity and the sense from ,,N” to ,,S” 

— the rotor rotated with œ, velocity by an external mechanical 
force. On the rotor winding turns are winded in which the 
electromotive force e, (EMF), is induced. The force can be 
described as follows: 


e= |k "B: sin(o,t) | (1) 


where: 

k — design coefficient of a given generator 

B — magnetic field intensity 

œ, — instantaneous angular velocity of generator rotor 


— the commutator, i.e. the ring fastened on the rotor and 
made of a conductive material. The ring is cut into 
segments which are electrically insulated from each other 
and form the so called commutator sectors (bars). To 
each of the sector the end of winding turn is connected; 
the commutator serves this way as a mechanical electric 
current rectifier 
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the electric brushes: „+” and ,,-”, which slide around 
commutator bars. To the brushes are connected wires 
conducting electric current to consumers. 


In order to increase magnitude of the inductivity B, rotor’s 


winding turns are placed on a core made of silicon steel sheets, 
that amplifies magnitude of EMF (e) about 10 thousand times. 
To correctly fasten the winding turns on the core they are placed 
in grooves. As cross-sections of the grooves have a toothlike 
form they are further called the ,, rotor teeth”. 


Comparison of the theoretical run described by Eq. 1 


(Fig. 2) with the real run of generator pulsation component 
(Fig. 3) does not show any similarity between them. 


Commutator 


Fig. 1. Rotor with two winding turns and commutator of four segments 


One rotation of rotor 


Fig. 2. Run of electromotive force between brushes in DC generator 


Groove pulsations of generator’s output voltage are 
produced as a result of change of reluctance due to whirling 
the grooved rotor. The groove pulsation frequency f, can be 
expressed as [1+4]: 


f,=Z-n/60 (2) 
where: 
Z — number of rotor grooves 
n — rotational speed. 


In the literature [1+3] voltage pole pulsations are often 
associated with the so-called rotational pulsations because of 
their mutual similarity. The phenenomenon of voltage pole and 
rotational pulsations is manifested in the run of output voltage of 
DC commutator generator, that can be observed in the form of 
changes of the run of the envelope shown in Fig. 3. Frequency of 
the modulation is directly proportional to the product of number 
of stator poles and angular velocity of rotor, whereas depth of 
its amplitude is proportional to changes of magnetic reluctance 
between rotor and stator. The pole pulsation frequency f, can 
be presented by means of the formula: 


f = 2p - n/60 6) 


where: 
p - 
Signal ofpole modulations carries information on anisotropy 


of sheets of generator magnetic circuit. In the subject-matter 
literature the pole modulation is usually associated with 


number of pairs of stator magnetic poles. 


Commutator Commutator sectors (bars) 
a) brush 
— 
Direction of motion 
Uk 
b) 
t 
c) 
| 
| 
Uz 
d) | 
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Fig. 4. Shape of pulsation curve for DC generator: a) developed view of 
mechanical elements of commutator unit, b) run of commutator pulsations 
- U; = f(t), developed view of rotor grooves, c) run of groove pulsations 
(continuous line) - U, = f(t) with added commutator pulsations (broken 
line), d) run of groove pulsation 


rotational modulation which is characteristic of that such 

modulation frequency is equal to the first or second harmonic 

frequency (in certain cases - to the first subharmonic) of 

rotational speed of generator rotor. The signal carries diagnostic 

information on the errors: 

— of workmanship of the generator, especially on inaccuracy 
of geometrical dimensions manifested as an asymmetry of 
air-gap between stator and rotor 


Angular (rotational) speed n = 8000 rpm/min; effective output voltage Uw = 28.53 V 
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Fig. 3. Changes of pulsation component of an aircraft DC generator under minimum load 
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— of assembling the generator, such as parallelism error, i.e. 
shift of rotor shaft axis with respect to that of drive shaft, 
sometimes called also eccentricity error, as well as angular 
error of shift of rotor shaft axis against drive shaft axis. 


Voltage commutator pulsations are associated with 
interaction of brushes and commutator. During armature 
rotation the brushes short-circuit alternately different number 
of winding turns, that introduces a change in number of turns in 
parallel branches and generates periodical pulsations of voltage 
at the brushes. The ferquency f, of the pulsations depends on 
number of commutator sectors and can be expressed by means 
of the formula [1, 2]: 


f.=K- n/60 (4) 


where: 
K — number of commutator sectors. 


GROOVE PULSATIONS 


The phenomenon of various groove pulsations is well 
described in the literature dealing with alternate-current (AC) 
induction generators [1, 2]. They do not possess any winded 
rotor and their useful signal is obtained from the stator winding. 
Their rotor is made of a ferromagnetic material (usually of 
a packet of silicon steel sheets) it has milled grooves (teeth) 
due to which modulation of magnetic field intensity of stator 
magnets is generated. In the generators groove pulsations are 
the crucial phenomena producing the useful signal. As there is 
no rotor winding - in contrast to the classical DC commutator 
generator - therefore only changeable component of pulsation 
is generated (commutator pulsations are not present because 
of lack of commutator and rotor winding). 

As results from the literature information [1, 2], to induction 
generators, in order to achieve an output voltage signal close to 
sinusoidal one, skew form of teeth is usually applied (Fig. 5b). 

Rotors having ,,dovetail” grooves (Fig. 5c) are rarely used 
in induction generators as then an unsymmetrical form of 
output voltage appears [1]. However such shape of grooves is 
commonly used in classical DC commutator generators. The 
shape makes firm fastening the winding onto rotor, possible. 
Time intervals between crossings of groove pulsations through 
the rotational speed reference level set for generator rotor, 
œ, = const., are dependent only on error in angle of milling the 
teeth. As the errors cyclically appear after every full rotation of 
the rotor they can be easily filtrated out. However the fact of stiff 
mutual angular position of grooves remains undeniable. Hence 
for œ, = var the time intervals between successive „zero”- level 
crossings (after filtrating any possible errors in milling the rotor 
grooves) will constitute a measure of instantaneous changes 
in angular velocity of rotor. The described features of groove 
pulsations have been used as a source of diagnostic information 
on technical state of generator drive system, on the basis of 
which the FDM-A diagnostic method (described in [5]) has 
been elaborated. 

The measuring of amplitude of groove pulsations makes it 
possible to localize breaks in rotor winding. As results from 
the data collected by these authors [7] after a winding failure 
a decrease of the relative value (related to the effective value 
of the generator output voltage U ) of the groove pulsation 6, , 
is observed. The value can be expessed as follows: 


=E (U aa Unnan) * 100%/Z-U, © 


max m 

where: 

m — natural number being that of successive interval of 
groove pulsation 
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Fig. 5. Typical runs of output voltage of induction generators with rotor 
teeth of the form: a) trapezoidal, b) rectangular, c) „ dovetail” — like 


— maximum instantaneous value of pulsation component 
within a given interval m 

— minimum instantaneous value of pulsation component 

within a given interval m 

number of rotor grooves. 


Simultaneously, after a failure of DC generator winding, 
the effective output voltage value changes, AU „(practically 
imperceptible, especially at lower values of œ, shown in 
Tab. 1 based on the data of [7], occur. The relative value of 
the changes, 6U, , due to failure of a single winding turn, 
does not exceed 0.01 %. In practice to detect a generator 
failure, i.e. occurrence of a break of its winding, by aircraft 
personnel under operation is entirely impossible. However 
this is fully possible, as results from experience gained by the 
team supervised by these authors, by making use of special 
measurement instruments. 

A failure of generator winding, e.g. its break, results in the 
decrease of the values of groove pulsation, given as the index 
Aò, in Tab. 1, from 0.8 % to 1.6 %, that can be practically 
measured by using measurement instruments of 0.1% class. 


POLE PULSATIONS 


The phenomenon of pole pulsations can be clearly observed 
on the run curve of DC generator output voltage [5] in the form 
of amplitude modulation, shown in Fig. 3. The modulation 
frequency is directly proportional to the product of number of 
stator poles and angular velocity of rotor, and the amplitude 
depth - proportional to magnetic reluctance changes between 
rotor and stator. The signal carries information on anisotropy 
of sheets of magnetic circuit of generator. The modulation 


Tab. 1. Parameters of groove pulsations before and after failure of generator winding 


State of winding 


Capable 


Failed 


Indices used 


to compare 
parameters 
before and after 


can introduce small errors in measuring AT,. It can be easily 
filtrated out because of its repeatability characteristic for 
a given generator. The relative value of pole pulsation, ò, , can 
be expressed as follows: 


+U 


ô, = iU U min 3 1 00%/U nax o min 6 } MAX (6) 


max 0 


natural number standing for successive number of 
pole pulsation interval 

maximum instantaneous value of voltage pulsation 
component in the o-th period 

minimum instantaneous value of voltage pulsation 
component in the o-th period. 


max o 


min o 


The pole pulsations carry a few kinds of diagnostic 

information: 

a) phase parameter informs on possible errors in geometrical 
distribution of stator pole shoes, 

b) pulsation amplitude (run of the envelope shown in Fig. 3) 
generally shows a non-uniformity of magnetic field 
distribution under stator magnetic poles and, in some cases, 
a shorting or break of rotor or stator winding: 

- if pole pulsation amplitude reaches, during the whole 
period, values uniformly increased and close to those 
of rotor groove pulsations, it means that one turn of 
winding is overloaded due to an increased leakance of 
its insulation or a partial fault to frame or between rotor 
winding turns in a given groove 

- ifpole pulsation amplitude reaches non-uniform values 
during the whole period, e.g. during one rotation of rotor, 
the peak value of the envelope undergoes a decrease; it 


failure 


means that one turn of winding is overloaded due to an 
increased leakance of its insulation or a partial fault to 
frame or between winding turns of one pole of stator 

- if pole pulsation amplitude value uniformly decreases 
during the whole period of rotor’s rotation, as shown 
in Tab. 2, it may constitute information on a break of 
rotor winding. 


With a view of a diagnostic complexity of such signal and 
of its small amplitude with respect to the carrier component 
(groove pulsation), location of failed winding turns by 
measuring pole pulsations seems rather inaccurate. 

However the pole pulsation signal becomes greatly 
increased in case of a failure, e.g. shorting of an arbitrary 
winding, consquently its amplitude increases many times with 
respect to that of groove pulsations. The shorting phenomenon 
is below described in detail as commutator pulsation amplitude 
greatly increases during shorting the winding turn. 


COMMUTATOR PULSATIONS 


The phenomenon of commutator pulsations has been not 
used in the FDM-A method [5, 6] as it has been regarded 
as a disturbing signal. The investigations performed under 
supervision of these authors [5, 8] have showed that amplitude 
value of the pulsations is directly proportional to current- 
load level. In Fig. 4 are presented mutual relations between 
commutator pulsations and groove ones, as well as location of 
rotor grooves and commutator bars. From the investigations 
performed with the use of an aircraft DC generator it results 
that at the generator’s current load lower than 10% of its rated 
value the commutator pulsation amplitude (Fig. 4b) is rather 


Tab. 2. Parameters of pole pulsations before and after failure of generator winding 


[rpm/min] | 4000 


9500 | State of winding 


Capable 


Broken 


Compari-son 
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unperceptible against groove pulsation background (Fig. 4d). 
At the current load of the order of 10% the pulsations are 
barely visible on the output voltage run. Angular displacements 
of particular halves of sinusoids of commutator pulsations 
(Fig. 4d) change with respect to groove pulsations and undergo 
individual angular displacements due to mechanical vibrations 
of brushes in brush-holder and during current loading the 
generator. Therefore the commutator pulsations cannot be 
used to diagnosing magnitude of failures of drive system’s 
kinematic pairs. 

Under rated load, the peak value of commutator pulsations 
reaches the level of about 50% of groove pulsations. It means 
that they may serve as a source of diagnostic information on 
e.g. commutator-brush unit failures. 

Trials of a controlled shorting in rotor have provided 
interesting data. In the case of shorting in the middle of one of 
the rotor winding turns it was revealed that the pole pulsation 
visible in Fig. 6 as the slow-varying component, became 
dominating, and the commutator pulsation visible in Fig. 6 as 
the fast-varying component, appeared to be that of the second 
order. 
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Fig. 6. Diagram of DC generator output voltage in the case 
of rotor windings shorting (the shorting localized 
in the middle of one of winding turns) 


However if the pole pulsation component was stable 
with respect to its frequency and amplitude, the commutator 
pulsation component reached its greatest value in the instant 
of passing under successive pole of generator stator. The 
unambigous change of the relation between amplitudes of 
pulsation components makes it possible to detect shortings in 
rotors of commutator generators. 


CONCLUSIONS 


— Inthis paper different kinds of output voltage pulsations of 
DC generator were described. Although they simultaneously 
occur in practice their amplitude-phase relations are very 
different. The subject-matter literature fully describes each 
of them separately. 
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— These authors, basing on their personal experience, have 
made an attempt to highlight practical relations between 
the pulsations. The pulsation component carries several 
diagnostic signals both concerning technical state of drive 
system and DC generator itself, i.e. the very source of the 
information, which has been not mentioned at all in the 
literature. 

— The diagnostic symptoms contained in the pulsations, 
precisely recognized by these authors, have been 
implemented by them to practice a few years ago. The 
other, ambiguous and not fully identified ones will be 
ready for application only after performing many arduous 
investigations aimed at finding accurate relations between 
successive parameters of kinematic faults and parameters 
of output voltage component. 
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ABSTRACT 


Article Deals with problem of naval fuels’ microbial contamination. This problem in Poland is not well 

known and described, but operational problems encountered very often prove that this is a very serious 

problem. Examples of laboratory tests performed on real samples are presented as well as possible effects 

of microbial contamination are discussed. Authors also suggest continuous monitoring of naval fuels as 
main preventive tool. Article covers following topics: 

e general information about microorganisms present in naval fuels 

e conditions necessary for microbial growth, with special attention to conditions encountered on vessels 

* operational problems resulting from microbial contamination 

e methods of microorganisms detection and examples of preventive actions. 


Keywords: naval fuels; microbes on fuels; microbial growth; microorganisms; fungi; ulphate-reducing bacteria 


INTRODUCTION 


It has been known for over 100 years that some 
microorganisms can feed with hydrocarbons but only for about 
40 years people have started to realise the adverse effect of 
microbes on fuels’ cleanness as well as their corrosiveness to 
materials used in distribution systems (storage and operational 
tanks, pipelines etc.) and in fuel systems of vessels. 

The core factor for microbial growth in fuels, not only the 
naval ones, is presence of water. Every case of fuel contaminated 
with water brings real and big risk of microbial growth. In fuel 
tanks ideal conditions for microbial growth exist on fuel-water 
surface. Fig. | presents, as a scheme, storage tank with layers 
that exist during fuel storage. 


Fuel 


Mix of fuel with biological layers 


Free water with particulate matters 


Sludge 
Fig. 1. Scheme of layers present in storage tank 


Microbial growth is stimulated also by other factors, i.e. fuel 
composition (as naval fuels most vulnerable are heavy fuels, 
fuel distillates (diesel fuels) and light fuels), storage conditions 
(temperature between 15 and 40 °C) and amount of free water 
in storage tank. 

Microbial contamination of naval fuels has posed some 
operational problems for many years. Such problems can be 
avoided by discarding of water, but practically it is impossible 
to discard all of the water. Water as separate phase in fuel tanks 
can origin from condensation of water dissolved in fuel as well 
as from water deliberately added to fuel tanks for ballast and 
to eliminate vapour space. 


MICROORGANISMS ENCOUNTERED 
IN FUELS 


Microorganisms having the biggest effect on fuels’ 
performace can be differentiated into two main categories: 
e fungi 
e sulphate-reducing bacteria. 


There are also yeast, but they are nor as important as fungi 
and bacteria. 

Fungi can create coherent mats at the fuel - water interface. 
Fungi enter the fuel from soil as airborne spores. 

Bacteria usually enter fuel tanks with water. They grow in 
water under the fuel layer, without oxygen presence [1]. 

Gloriah Hettige in her work [2] has isolated thirty one types 
of fungi and five types of bacteria. Neyhof and May have also 
described other types of fungi and bacteria [3]. 

The main types of fungi were the following: 
e Cladosporium resinae 
e Penicillium corylophilum 
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e Paecilomyces variotti 
e Aspergillus 

e Mucorales 

e other, unidentified. 


The bacteria were the following: 
e Pseudomonas 


e Desulfovibrio. 


Among yeast there were the following types: 


e Candida 
e Rhodotorula. 
MICROBIOLOGY OF FUELS 


The first reports regarding microorganisms were in the early 
1930’s. The reports described problems connected with fuels 
and oils infected by bacteria [4, 5]. 

Most of work concerned aviation fuels, though probably 
there also were problems with naval fuels, but they were used 
in less sophisticated (compared to jet ones) engines. 

Liggett [6] has performed comprehensive work on microbial 
contamination of marine, rail and road diesel fuels. 

Hydrocarbons C,, - C,, are much readily assimilated by 
microorganisms than the C, - C, ones. Test results showed that 
microbial growth had taken place also in naval fuels, but effect 
of their activity was no as severe as for aviation fuels. 

This was why some contamination and corrosion in naval 
fuel recognised as of no importance so it was not justified to 
take appropriate countermeasures similar to ones for aviation 
(water discarding, microbial monitoring). 

As a result the problems regarding microbial growth at 
vessels became bigger. 

Such unfavourable occurrences were linked mainly with 
chemical changes in fuel resulting from refinery processes [7] 
as well as from wider use of additives, which sometimes 
promote microbial growth. 

Hill [8] has found severe microbial growth both in heavy 
fuel oils and lighter diesel oils. 

The biggest problems were because of fungi “Cladosporium 
resinae ” which were superseded by a lot of other types of fungi, 
bacteria and yeast. 

Distillate fuels contain ample of oxygen sustaining aerobes 
growth, but in case of dormant tanks, content of oxygen 
decreases and anaerobes such as sulphate-reducing bacteria 
Desulfovibrio can grow. Such bacteria never exist alone, but 
also together with other organisms, creating big complexes 
living in water and feeding on fuel. Fuel with its additives 
is not the only nutrition source. It is permanently used and 
replenished, opposite to water and microorganisms contained 
there. 

Conditions in naval fuel tanks, especially presence of large 
amounts of water, effect on variety of types of microorganisms 
as well as on their growth rate. It can be assumed that microbial 
growth in vessel’s fuel systems is bigger and faster than in 
aviation ones. 

The ratio of water to fuel in water-displaced tanks can 
be high. There is wide range of nutrients in such mixture, in 
addition to those already contained in water. The nature and 
quantity of such nutrients depend on type of water used to fill 
the tanks. It can be saline deep-seawater as well as estuarine 
one. Extent of colonisation can be affected by sewage content 
and trace metal contamination in the water. 

The ratio of water to fuel in standard fuel tanks is low 
and there is less particulate organic matter when compared to 
water-displaced tanks. Amount of nutrients will be determined 
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by the previous history of the fuel and such factors as condition 
of tank linings and seals. 


RESULTS OF MICROBIAL ACTIVITY 


Inconsiderable amount of water, even strongly infected, 
does not influence on significant fuel deterioration. Such water 
can result in severe corrosion. It is well known that more 
the water much bigger threat to cleanness of fuel in terms of 
microbiology. Activity of microorganisms can lead to: 

e presence of slimes and clogging pipes, valves and filters 

e coalescers malfunction due to fungi growth on filter cloth 
sleeve 

e malfunction of volume measurement systems 

e accelerated corrosion due to aggressive substances (acids, 
sulphides) and electro-chemical corrosion cells creating, 
corrosion inhibitors depletion and degradation, protective 
coatings degradation and preventing formation of stable 
oxide layers. 


Symptoms of severe contamination of naval fuels are the 
following: 
e filter clogging 
e coalescers malfunction 
e injector fouling 
e fuel system components corrosion. 


MICROORGANISMS DETECTION 


Crucial matter is to obtain water bottom sample. For such 
sample it is possible to perform both qualitative and quantitative 
testing using a routine culture methods. But sometimes it 
is impossible to obtain such excellent sample, so in order 
to evaluate level of fuel contamination we can use different 
methods. Some of them are summarised below: 


Microscopic method 


It is used mostly to detect motile bacteria, fungal hyphae, 
protozoa, algae and unusual debris. It also allows identification 
of small particles such as non-motile bacteria and spores in 
the presence of large amounts of mineral particulate matter. 
Disadvantage of this method is the high cost of devices and 
necessity of highly qualified personnel, but major advantage 
is short time of analysis. 


Centrifuge method 


It bases on centrifuging sludge sample previously 
homogenised. The length of specific layers in test tube is read 
from calibration graph. Its disadvantage is, as in previous 
metod, high cost of equipment. 


Chemical methods 


e pH of water phase measurements 
e acidity measurements 

e metallic elements content 

e sulphides presence 

e burning at 500°C. 


Membrane filters method 
It is based on passing fuel being examined through defined 


membrane filter and evaluation of substances collected on 
filter’s surface. 


Dip-slide method 


The nature of this method is to dip a small plastic paddle 
coated with nutrient into a sample. Next step is to return the 
paddle into its container and to incubate. One type of nutrient 
detects ant evaluates bacteria while the another detects yeast and 
fungi, and allows for quantitative assessment of contamination. 
Disadvantage of this method is long incubation time, up to 
3 days, but major advantage is simplicity so it can be used my 
not qualified personnel. Examples of results obtained with this 
method are presented on Fig. 2. 


Fig. 2. Plastic paddle coated with nutrient (prod. Merck) 
used in dip-slide method 


ATP method 


This is a new method evaluated recently. It uses presence 
of ATP in microorganisms contained in fuels. Its main 
advantage, when compared with other methods, is short time 
of measurements that allows for reliable results within approx. 
10 min. Example of test device used in this method is presented 
on Fig. 3. 


Fig. 3. Hy-Lite device (prod. Merck) used in ATP metod for microbial 
contamination assessment. 


EXAMPLES OF NAVAL FUELS 
MICROBIAL CONTAMINATION 


Examples of naval fuel samples, where microbial 
contamination has been confirmed, are presented on Fig. 4. The 
samples were taken from operating vessels. Easily visible is fuel 
colour change and microbial layer on fuel-water interface. 

Slides with cultures of fungi, bacteria and yeasts are 
presented on Fig. 5. 


Fig. 4. Naval fuels samples with confirmed microbial contamination 


Fig. 5. Slides after tests of microbial contamination 
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COUNTER MEASURES AND 
PREVENTIVE ACTIONS 


Fuel composition 


Fuel should be formulated that it is no longer good nutrient 
for microbial growth. Straight hydrocarbon chains are more 
vulnerable than branched ones. Some additives promotes 
microbial growth and some prevents that. So far, there is 
no comprehensive study concerning fuel composition effect 
on microbial growth, and refineries are using established 
production processes, therefore changes in fuel composition 
are not currently feasible. 


Good housekeeping 


This mean includes: 

e minimising water content 

e rust and scale discarding 

e high temperatures avoidance 

e avoidance of fuel contact with contaminated surfaces 

e active fight against microbial growth at the beginning of 
the process of growing. 


Mechanical decontamination 
This preventive action includes: 
e filtration 
e centrifugation. 
Biocide treatment 
Biocides are the chemicals used to destroying microbes. 


There are wide range of chemicals used as biocides. Some of 
them are the following: 


¢ Biobor JF 

e Kathon 886 
e Grotan OX 
e Bodoxin 


e Omadine TBAO 
e Bioban FP etc. 
e Grota Mar. 


Biocides should have the following properties: 

e be combustible, without ash 

e have no adverse effect on fuel properties or engine 
performance 
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e be soluble or very miscible in fuel, but preferably in 
water 

e be not corrosive 

e be able to total killing of all microbes at minimum dose. 


CONCLUSIONS 


e Microbial contamination of naval fuels pose problems and 
threat to proper use of vessels. 


e Effect of microorganisms can be minimised by adhering to 
rules of good housekeeping, as well as by the use of relevant 
chemicals (biocides and biostats). Use of the chemicals 
should be based on relevant laboratory testing to check its 
influence on fuel system’s components. 
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ABSTRACT 


This paper presents a description of new global monitoring system for containers with its layer-modular 

structure, as a solution for enhance security and efficiency of container transport with particular emphasis 

on the practical implementation of that system for maritime container terminals. Especially the Smart 

Container Module (SCM) architecture and its operation as a part of the Self-Organizing Container 
Monitoring Network is presented. 


Keywords: wireless monitoring; container transport; container module realization; 
global system concept; ad-hoc networks 


INTRODUCTION 


Since 26 March 1956 when the first goods transport was 
made by using freight containers on board of the Ideal X 
— the world’s first container ship - this kind of transportation 
gained very rapidly on popularity. Meanwhile, the concept of 
the container has been clearly defined, its dimensions have 
been standardized and on the whole world have been created 
numerous container terminals: land and sea, specialized in 
handling this type of cargo. Currently container transport 
includes more than 90% of the global trade and takes place 
along trade lanes, both marine and terrestrial. 

However, despite all the advantages of this mode of 
transport, it carries a lot of risk, as a result of the exploitation 
of containers. First of all, containers are closed during the 
transportation, so the load on their interior remains beyond 
any control. For this reason, in peer-reviewed literature on 
the safety of transport, container is known as the Trojan 
Horse of twenty-first century. In this state of affairs it is 
particularly important issue of controlling and monitoring 
of cargo container, both during transportation and storage. 
Several years ago, the Group of Self-Organizing Ad-Hoc 
Wireless Sensor Networks was founded at the Department 
of Radiocommunication Systems and Networks in the 
Gdansk University of Technology. This Group developed the 
original concept of a global wireless system for monitoring 
of container cargo, with particular emphasis on the practical 
implementation of that system for maritime container 
terminals [1, 2, 3]. 


A GLOBAL GRASP OF THE SYSTEM 


The proposed new monitoring system solution for containers 
has a layer-modular structure (Fig. 1). The basic layer of the 
system is the Smart Container Module (SCM), designed for 
installation on containers and eventually should be an integral 
part of each container. Another, higher layer of the system is 
the Ship/Port Subsystem with the self-organizing container 
monitoring network, based on the 802.// standard, working 
on the container ship or in the container terminal and the sea 
ports. Due to the fact that the typical container ship can carry 
up to dozen thousand containers, there is a necessity for the 
location a special database on the ship. In case of container 
terminals and sea ports there is the Port Administration Data 
Base, equipped with wireless access points located at entrances 
to ports and wireless radio networks covering container depot. 
All containers data are regularly sent to the ship/port network 
controller. Information can also be sent on request. All of the 
Port Subsystems located inside the country are connected to 
each other and provide all data to the Land Subsystem - The 
State Administration Data Base. The Land Subsystem is the 
third layer of the monitoring system. 

Transmission of information between each part of the 
system is in accordance with certain rules, using an appropriate 
transmission medium. In the case of a container ship, the period 
of reporting depends on its position, i.e. on distance to the 
mainland, the home port and/or destination port. Frequency 
of sending the relevant data increases with the approach to 
the mainland. 
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Fig. 1. Architecture of the global wireless monitoring system for containers 


The Global Subsystem is the highest layer of the monitoring 
system, which can use the database of a new LRIT (Long 
Range Identification and Tracking) system for long-range ship 
monitoring. Data are transferred to the global database via 
satellite system (e.g. INMARSAT), when the vessel is on the 
open sea. Approaching the mainland and entering the harbour, 
information are communicated to the Port Subsystem that 
supports appropriate area. Transmission of this information 
is carried on using the GPRS system. It is also necessary to 
realize co-operation and data exchange between local and 
global databases, and also co-operation between databases 
inside the country. It can be realized by means of fixed lines, 
e.g. optical fiber. 

In case of land transportation all data about monitoring 
cargo are sending via GPRS to State Administration Database. 
The SCM modules are mounted on containers carried on 
container lorries and there is no need to organize SCM modules 
in ad-hoc network. 

The recipients (users) form a separate layer of the system. 
We can distinguish two classes of users: the services responsible 
for cargo security and the carriers. The first group is formed 
by the relevant services, such as Coast Guard, Customs and 
Police. The owners and the shipping companies belong to 
second group of users, which access to the system is carried 
out by authorized employees. 


THE SMART CONTAINER MODULE 


The Smart Container Module is one of the most important 
elements of the monitoring system [1]. Fig. 2 depicts its 
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functional diagram. The versatility of the SCM stems from the 
possibility of parameters measuring and monitoring of state 
within the container. To this end, the module is equipped with 
various sensors placed inside the container [4]. This sensor 
network enables detection of movement, temperature changes, 
pressure, humidity, presence of gases or radioactive substances, 
etc. An important aspect is also the location of the container, 
so each SCM is equipped with the GPS receiver. The GPRS 
module and the WiFi module are components of SCM. This 
parts of container module enable communication with port 
database (GPRS) or the local wireless network working on the 
container ship (WiFi). 

Transportation documents concerns shipment (such as 
contents of the container, a shipper, a point of shipment, 
a point of destination and other information) are stored in 
SCM’s memory. This data are sent to the appropriate database, 
whenever necessary. In case of lack of connection with the radio 
resources of the system, it is necessary to store data about status 
of container in the module’s internal memory and send them 
shortly after the SCM gets in a range of network. 

The container module has to be characterized by very 
low power consumption to enable unattended work for 
a long time (from several months to even few years). 
Equally important is compact construction of the SCM and 
to minimize its dimensions. The module’s casing has to be 
matched to construction of the container, to damage prevent 
during transshipment. It also has to be resistant to the weather 
conditions and a protection against removal or modification 
of its interior, with automatic alert generating in these cases 
is needed. 
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Fig. 2. Functional diagram of the Smart Container Module 


During the research on the SCM, it was decided to use two 
types of industrial computers: the Axiomtek SBC84710 and the 
Advantech ARK1380. Three prototypes of the SCM have been 
produced basing on these computers. The Smart Container 
Module prototype using the ARK1380 is shown in Fig. 3. The 
module is enclosed in a sealed casing and there are various 
connectors available, as following: 

e COMI, COM2 - used to connect Sensors’ Controllers in 
various configurations 

e USB], USB2, USB3, USB4 - two of them are used to connect 
external GPRS and GPS modules, next two are designed 
for future purposes 

e LAN - used to connect the SCM with the ship/port local 
network, in case of usage the SCM as a gateway 

e DC IN - used to connect the supply voltage range from 
9V to 35V, but because of the usage of this voltage by 

Sensors’ Controllers, it is recommended not to exceed 

a value of 15V 
e PCMCIA - used to connect the WiFi module 
e VGA - used to connect the monitor 
e PS2 - used to connect the keyboard 
e AUDIO - not used, designed for future purposes 
e LVDS - not used, designed for future purposes. 


Additionally, there are LEDs placed on the SCM’s casing, 
which indicate supply voltage connection and work of the hard 
disk. There is also the on/off button placed on the casing. 

The internal construction of the Smart Container Module 
is presented on an example of the SCM prototype based on 
SBC84710. The top view of internal construction of the SCM, 
with marked available interfaces, connectors, power supply 
voltage and connected the WiFi module is shown in Fig. 4. The 
bottom view of the internal SCM construction, with marked 
the RAM chip and the Compact Flash Card, which stores 
the operating system and all of the SCM’s data is shown in 
Fig. 5. 

The container module is powered by two batteries thus it 
enables unattended work for a long time, which depends on 
batteries capacity. The first of batteries has a smaller capacity 
and it is enclosed in the SCM’s casing. The second of batteries 
could have as large capacity as necessary. It is placed inside the 
container and connected to SCM. The SCM is equipped with 
a voltage control circuit, whereby it is possible to distinguish 
SCM turning off for lack of power against other cases. 

The Sensors’ Controller have been designed in a manner 
making construction of the measuring module independent from 
a type of using computer model and to enable free deployment 


Fig. 3. A prototype of the Smart Container Module 
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Fig. 5. Bottom view of the internal SCM construction 


of reconfiguring sensors inside the container. The exemplary 
Sensors’ Controller board is shown in Fig. 6. Composition of 
this circuit includes three parts: the microcontroller, the RS485 
converter with auxiliary circuits and the set of sensors dedicated 
to the specific cargo container. Sensors, which can be connected 
the Sensors’ Controller with: 


Fig. 6. The sensors’ controller board 
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e sabotage - this sensor allows to detect a sabotage, e.g. 
SCM’s casing tamper or cutting off any of Sensors’ 
Controllers 

e door switch - this sensor indicates whether the door of the 
container is open or not 

e motion - this sensor allows to detect movement inside the 
container 

e smoke - this sensor allows to detect smoke inside the 
container 

e temperature - this sensor measures ambient temperature 
inside the cargo container 

e humidity - this sensor measures relative humidity inside the 
cargo container 

e acceleration - this sensor measures acceleration on three 
axes (range to 10g or to 50g) 

e any other sensor with electrical output. 


Usage of Sensors’ Controllers allows connect many 
sensors to the SCM simultaneously. Additionally it allows to 
unrestricted changes types of sensors and their deployment 
inside the container. 


SELF-ORGANIZING NETWORK 
OF SCM MODULES 


Because of the vertical way of containers storage on 
container ships, there is a problem of access to the on-board 
wireless network for Smart Container Modules mounted on 
the containers at the lowest level. Additionally, GPS satellites 
are beyond the sight of view of the GPS receiver, installed in 
the adversely located SCM modules, thus it is impossible to 
read the geographical position. In order to solve this problem 
SCM modules have been programmed in a manner that allows 
them to self-organize in ad-hoc network. During the research 
a new algorithm for data transmission control over a multi- 
hop ad-hoc network (assuming slow-moving nodes) has been 
developed [3]. 

By dint of this solution it is possible to transfer data about 
cargo to the Ship Subsystem’s database, even from the Smart 
Container Module installed on the lowest container. The data 
about the cargo are transferred to adjacent modules by the WiFi 
radio link. If module has a connection to the on-board wireless 
network, it will forward data to the database. If this module 
is beyond the reach of the on-board wireless network, it will 
forward the data to the next SCM, etc. Additionally, the self- 
organizing ad-hoc network enables SCM modules beyond the 


Fig. 7. The SCM self-organizing ad-hoc network tests 


reach of GPS system to read and save the geographical position 
on the basis of data from modules at the top of containers’ stack. 
The ability to self-organize modules in the ad-hoc network can 
also be useful in harbours or container terminals, where direct 
access to the Port Subsystem’s wireless network or receiving 
signals from GPS satellites can be impossible. Fig. 7 depicts the 
laboratory for testing the Self-Organizing Container Monitoring 
Network composed of three Smart Container Modules. 

The results of these tests have proved the ability of testing 
network to self-organize in real time, especially to choose SCM 
module to communicate with the database. 


USER INTERFACE 


In the course of the research, the user interface including 
a database storing all data about the system has been created. 
Depending on the granted access privileges, user may get 
the following information: state of the system, defined 


types of containers, registered containers (with the assigned 
sensor network), ongoing and completed shipments and 
communication history. In addition, there is the ability to users 
manage (adding/removing user, changing password) only for 
users with administrator privileges. 

Each transport can be remotely monitored. Cargo 
information (parameters from the sensor network, current 
position and route of transport) are available in three forms: 
a table, graphs and a map. Fig. 8 depicts a map of portion of the 
transport route. This function uses the Google Maps available 
for free by Google. 

Data can also be presented in a transparent form of graphs, 
which are created on the basis of data from various sensors. 
Sample record of the temperature inside the cargo is presented 
in Fig. 9. 

Using the interface is intuitive and user-friendly, which 
greatly facilitates access to the system at any point on earth 
with access to the Internet. 
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Fig. 9. Sample view of the user interface in the control of temperature scope 
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CONCLUSIONS 


This is a fully operational global system for wireless 
monitoring of cargo containers that enhances the safety 
and effectiveness of this mode of transport. The Smart 
Container Module (SCM) is a basic element in a global 
grasp of the system. The module is equipped with various 
sensors placed inside the container that enable detection 
of many parameters. Each SCM is equipped with the GPS 
receiver, the GPRS module and the WiFi module. Due to 
the vertical way of containers storage, there are problems 
with individual SCM to communicate with the local 
wireless on-board network or in the container terminal. 
The self-organizing network of smart container modules 
is a proposed solution of those problems. 

Described concept of global system for monitoring of cargo 
containers, especially self-organizing monitoring network 
consisting SCM modules is a part of a wider security 
subject matter, which currently has a high priority, both in 
Poland and the world. In addition, increasing the efficiency 
of container transport is possible by dint of the safe and 
efficient user interface. 
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ABSTRACT 


This study provides a framework for integrating various planning activities in container terminals. First, 

we introduce various planning activities in container terminals and identify decision-making problems 

for each planning activity. Input parameters, decision variables, objectives, constraints, time buckets, and 

the planning horizon for each decision activity are identified. Next, we introduce the concept of a resource 

profile and planning procedure simultaneously by considering availabilities of various resources and 
resource requirements in a planning activity. 


Keywords: container terminals; planning; resource capacity 


INTRODUCTION 


Port container terminals are located at places where 
containers are transshipped from a transportation mode to 
another. Their main functions are to provide transfer facilities 
for containers between vessels and land transportation modes, 
such as trucks and trains. They are characterized by highly 
complex systems that involve numerous pieces of equipments, 
operations, and handling steps. Operations in container 
terminals can be classified into a vessel operation process 
during which containers are discharged from and loaded onto 
a vessel and receiving and delivery operation processes during 
which containers are transferred from and to external trucks. 
During these operations, assigning resources to these operations 
and scheduling these operations become major planning issues 
in container terminals. 

Many researchers reviewed planning problems in an 
operation of container terminals (Ramani, 1996; Bontempi et 
al., 1997; Meersmans and Dekker, 2001; Vis and de Koster, 
2003; Steenken et al., 2004; Murty et al., 2005; Crainic and 
Kim, 2007). Ramani (1996) divided the basic task in the 
management of container terminals into a berth allocation, 
yard planning, stowage planning, and logistics planning in 
container operations. A berth allocation issue is to plan which 
berth is to be assigned to a given ship for loading and unloading 
its containers. A yard planning involves the allocation of 
storage spaces to import, export, and transshipment containers. 
A stowage planning assigns stowage locations to outbound 
containers in bays of a ship. A logistic planning deals with 
scheduling and coordinating the operation of port equipments, 
such as quay cranes, prime movers, and yard cranes for 


moving containers among different sources and destinations 
(for example, gates, vessels, rail stations, storage yards, and 
container freight stations). 

Bontempi et al. (1997) assigned a different time horizon 
to each planning problem in container operations. They used 
the solution obtained in long-term problems (container storage 
policies) as an input for mid-term problems (resource allocation 
problems) and the solution obtained in mid-term problems 
(resource allocation policies) as an input for short-term 
problems (load and unload scheduling problems). Meersmans 
and Dekker (2001) and Vis and de Koster (2003) distinguished 
decisions on container handling operations into strategic, 
tactical, and operational levels according to the time horizon 
involved. A time horizon in decisions for the strategic, tactical, 
and operational level covers one to several years, a day to 
months, and a day, respectively. 

Steenken et al. (2004) reviewed terminal logistics and 
optimization methods. They described the important processes 
in container terminals that can be optimized by means 
of operations research methods: ship planning processes 
(consisting of berth allocation, stowage planning, and crane 
split), storage and stacking logistics, transport optimization, 
and simulation systems. Murty et al. (2005) introduced nine 
decisions to be made in daily operations: allocation of berths 
to arriving vessels, allocation of quay cranes to docked vessels, 
appointment of arrival times to external trucks, routing of 
trucks, dispatch policy for trucks at terminal gatehouses and 
docks, storage space assignment, yard crane deployment, 
internal truck allocation to quay crane, and optimal internal 
truck hiring plans. Giinther and Kim (2006) divided planning 
and control levels in container terminals into three categories: 
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terminal design, operative planning, and real-time control. 
A terminal design level contains multi-modal interfaces, 
terminal layout, equipment selection, berthing capacity, and 
IT-system and control software. An operative planning level 
contains storage and stacking policies, crane assignment and 
split, berth allocation, and stowage planning. A real-time 
control level contains landside transport, quayside transport, 
slot assignment, and crane scheduling operation sequencing. 
Crainic and Kim (2007) introduced models for the operational 
planning control in container terminals. Their models include 
scheduling of berths, scheduling of quay cranes, stowage 
planning and sequencing, storage activities in a yard, and 
allocation and dispatching of yard cranes and transporters. 

This study provides a framework for integrating various 
planning activities in container terminals. First, we introduce 
planning problems in container terminals and identify a decision 
activity for each planning problem. Input parameters, decision 
variables, time buckets, and the planning horizon for each 
decision activity are identified. Next, we introduce the concept 
ofa resource profile and planning procedure simultaneously by 
considering resource capacities and resource requirements. 

Section 2 provides a framework for the operation planning 
process in container terminals. Section 3 discusses resource 
profiles for certain decision activities related to various 
operational plans. Three subsections illustrate examples of 
the capacity planning in berth planning, scheduling of quay 
cranes, and yard planning, respectively. Finally, Section 4 
presents a conclusion. 


FRAMEWORK FOR A PLANNING 
PROCEDURE 


Container terminals perform various handling operations by 
utilizing such resources as quays or berths, quay cranes (QCs), 
storage yards (SYs), yard cranes (YCs), traveling areas (TAs), 
and transporters (TRs). A traveling area represents a traffic zone 
or a set of transfers for trucks. Two most important performance 
measures in a container terminal are turnaround time of a vessel 
and road trucks in a terminal. Turnaround time of a vessel and 
road truck highly depends on the capacity of applied resources 
and methods to allocate the capacity of resources in handling 
tasks. This section describes how the capacity of resources 
can be explicitly considered during each of planning processes 
and information on the capacity can be shared among different 
planning processes in container terminals. 

Resource capacities are represented as follows: the capacity 
of berths can be represented as the product of one dimensional 
space (usually length) and time. The capacity of QCs, TRs, 
and YCs are measured in QC times, TR times, and YC times, 
respectively. If all QCs have the same capacity, we can evaluate 
the capacity of QCs by multiplying the number of QCs by their 
available time. In a similar way, the capacity of other resources 
can be evaluated. Tab. 1 summarizes units to represent the 
capacity of each different type of resources. 


Tab. 1. Resources and units 


Berth (H) Length of the berth x berthing duration 
QC (C) 


Number of QCs x operation time 


TR (R) 
YC (Y) 
TA (A) |Number of TRs passing the TA x operation time 
SY (S) 


Number of TRs x operation time 


Number of YCs x operation time 


Number of slots (in TEU) x storage duration 
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All planning activities for operations must check the 
availability of resources required for the operations. An 
amount of resources required for an operation can be 
estimated as follows: berthing a vessel requires the resource 
of berths based on the length of a vessel multiplied by the 
occupation time. Unloading or loading containers consume 
the resource of QCs by the number of containers multiplied 
by the standard handling time per container. The workload 
in unloading and loading operations on TRs can be evaluated 
by multiplying the number of containers with the average 
transportation time per container including empty travels. 
The workload of unloading and loading operations on YCs 
can be calculated by multiplying the number of containers 
with the standard handling time of a YC per container. The 
workload on TAs represents the expected future occupation 
of TAs by TRs. The storage space of a SY requires a certain 
amount of reservation before the storage of containers and 
actual occupation by containers. Fig. 1 shows a planning 
procedure in container terminals. We classified this procedure 
as long-term, mid-term, and short-term according to the 
planning horizon. 


Long-term calling 
schedule 


= 


Long-term Berth planning 


Quay crane 
allocation 


a 


Quay crane split 


| 


Unload / load 
scheduling 


Yard planning Operator scheduling 


Mid-term 


Remarshal planning 


Real-time operation 
scheduling for QCs, 
YCs, and TRs 


Storage location 
assignmanet 


Short-term 


Fig. 1. Framework of the planning process 


Every planning process must consider the availability 
of related resources. Fig. 2 shows several key planning 
processes and their related resources that must be checked 
before commitment. The berth planning is a decision process 
on the berthing location and time for ships. A certain zone 
in a berth is assigned to a ship for loading and unloading 
containers for a certain period of time. Scheduling of QC 
works (split) is a decision process on a service sequence of 
bays in a ship by each QC and time schedule for services. 
Several QCs are usually assigned to one ship. Yard planning 
is a decision process on storage locations for containers during 
unloading operation, receiving operation, and remarshalling 
operation. Scheduling of unloading is a decision process on 
a work sequence and time for inbound containers that are 
discharged from a ship to storage locations in a SY. Scheduling 
of loading is a decision process on a work sequence and time 
for outbound containers that are loaded from a SY to storage 
locations in a ship. Remarshal planning is a decision process 
on movements of containers from a storage block to another 
and time of movements. 


Berth 
capacity 


Truck 
capacity 


Traffic area 
capacity 


Yard planning 


| Yard crane |<" 
capacity 


Remarshal 
planning 


| Storage yard ¥ 
capacity 


Fig. 2. Various operational plans and their related resources 


Berth planning 


Decisions to be made by each planning process are 
summarized in Tab. 2. The first line of each planning process 
represents the activity or activity unit on which a decision will 
be made. The second line represents the reference moment 
of the activity from which the time-phased consumption of 
resources resulting from the decision will be estimated. The 
third line represents contents of the decision to be made by 
each planning process. 

The following sections will discuss the resource profile in 
major decision processes for long-term and mid-term plans 
in container terminals: berth planning, scheduling of QC 
works, and yard planning. This study considers six resources: 
berths, QCs, TRs, YCs, TAs, and SY. The operation time 
applied in resources is configured as berthing time, handling 
time by QCs, transportation time by TRs, handling time 
by YCs, occupation time of TAs, and storage time of SY. 
For describing the resource profile and plan, the following 
notations are used: 


Quay crane 
capacity 


QC work 
scheduling (split) 


Unload / load 
scheduling 


Tab. 2. Contents of the decision for various operational plans 


Long-term 


Activity to be planned 


Berth planning (B) 


Berthing of each vessel 


Reference moment of the activity 


Beginning of the berthing 


Contents of the decision 


Berthing position and time of each vessel 


QC work scheduling (Q) 


Activity to be planned 


Loading or unloading task on deck or in hold of a bay by a QC 


Reference moment of the activity 


Beginning of each task 


Contents of the decision 


Schedule for QCs to discharge (or load) containers from (or onto) vessels 


Yard planning (P) 


Activity to be planned 


Receiving outbound containers for a vessel or unloading inbound containers by a QC 
for a vessel for a period 


Reference moment of the activity 


Starting of arrivals of outbound containers at a gate or unloading of inbound containers 
from a vessel 


Contents of the decision 


Storage positions for outbound containers for a vessel arriving during a period or 
inbound containers discharged from a vessel 


Remarshal planning (M) 


Activity to be planned 


Moving a set of containers from a block to another for a period 


Reference moment of the activity 


Starting of movements 


Contents of the decision 


Containers to be moved and their sources and destination positions for a period 


Short-term 


Activity to be planned 


Unload scheduling (U) 


Unloading containers on deck or in hold of a bay 


Reference moment of the activity 


Starting of unloading 


Contents of the decision 


Discharging sequence of individual inbound containers 


Load scheduling (L) 


Activity to be planned 


Loading containers onto deck or into hold of a bay 


Reference moment of the activity 


Starting of loading 


Contents of the decision 


Loading sequence of individual outbound containers 
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Indices 

r — index used in resources where r= H (berth), C (QC), 
R (TR), Y (YC), A (TA), and S (SY) 

t — index used in periods where t = 1, 2, ..., m 

a — index used in activities where a = 1, 2, ..., n. 

Problem data 

si, — aunit amount of resource r that must be used with 
the time offset of t for carrying out activity a. In 
a berth planning, for example, if a vessel is decided 
to berth at a quay at period p (let this be activity “B’’), 
the operation time of QCs at period (p + t) will be 
required by the amount of sho. 

Sets 

T, —  asetoftime offsets in which resource r is consumed 
by activity a 

B — a set of activities related to berth plans. Each 
activity corresponds to a decision on the berthing 
of a vessel 

Q — aset of activities related to QC work schedules. Each 
activity corresponds to a decision on the work for 
a vessel by a QC 

P — aset of activities related to yard plans. 


Decision variables 

X, — a decision vector for activity a. As an illustration, 
a decision vector in a yard plan for outbound 
containers is a set of storage blocks for the arriving 
outbound containers to be stacked and the amount 
of containers to be stored at each storage block. 


BERTH PLANNING 


A berth planning determines the berthing position and 
time of a vessel. Tab. 3 shows input parameters, decision 
variables, objectives, and constraints for the berth planning. 
A berth planning requires such input data as a calling schedule, 
favorable berthing location, length of a vessel, required draft 
of each vessel, required unloading and loading time for each 
vessel, and various constraints and priorities for locating and 
sequencing a vessel. One of the objective in a berth planning is 
to locate vessels at the most favorable position on the quay that 
will reduce the container delivery time between the marshaling 
yard and QCs and also to make vessels to depart the port before 
their committed due times. Also, the waiting time of vessels at 
a port must be minimized. 


Tab. 3. Definition of problems in a berth planning 


Calling schedule of vessels 

Favorable berthing location of each vessel 
Length of each vessel 

Draft required for each vessel 

Number of unloading and loading containers of 


Input 
parameters 


Decision |Berthing position of each vessel 


variables |Berth time of each vessel 
To minimize delays in the departure of a vessel 
To minimize the travel distance between the 
shore and the yard for all containers in a vessel 
To minimize the waiting time of a vessel at a port 
Depth of water for berths 
Constraints|Arrival time of a vessel at a port 

Availability of resources 


Objectives 
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Resource requirements depend on the number of unloading 
and loading containers and are evaluated by using the load 
profile of each resource in a berth planning. The requirement 
of each resource can be determined as a time-phased way with 
respect to the berth time of a vessel. Decision variables in 
a berth planning consist of the berthing position and berthing 
time of each vessel. By adding the schedule of a vessel in 
a berth plan, various resources are required as shown in 
Figs 3-(a) and (b). Figs 3-(a) and (b) show the resource profile 
for outbound and inbound flows, respectively. 

A resource profile can be evaluated by using several 
different ways: s3,, indicates the amount of berths required 
in the berth plan B at a day when a vessel arrives at a berth. 
It can be evaluated by using (the length of a vessel plus the 
allowance between adjacent vessels) x (the expected berthing 
duration of a vessel). s,, can be evaluated by using {(the time 
for a QC to transfer an outbound container to a slot of a vessel) 
x (the percentage of loading containers among all containers for 
a vessel) + (the time for a QC to transfer an inbound container 
to a TR) x (the percentage of unloading containers among all 
containers for the vessel)}. Standard time that includes not only 
pure operation time but also unavoidable delays and personal 
and fatigue allowances must be used. 

Because TRs are required for loading and unloading 
operations, Ss}, can be evaluated in the similar way as s}\.. Shy 
for t <0 can be evaluated by using (the time for a YC to receive 
an outbound container from an external road truck) x (the 
percentage of containers, among all the outbound containers, 
arriving at a SY on the t day before the berthing of a vessel). 
Shy for t > 0 can be evaluated by using (the time for a YC to 
transfer an inbound container to an external road truck) x (the 
percentage of containers, among all the inbound containers, 
leaving a terminal on the t" day after the departure ofa vessel). 
Shy for t = 0 can be evaluated by using {(the time for a YC to 
transfer an outbound container to a TR) x (the percentage of 
loading containers among all containers for a vessel) + (the 
time for a YC to receive an inbound container from a TR) x 
(the percentage of unloading containers among all containers 
for a vessel)}. s$, for t < 0 can be evaluated by using (the 
cumulative percentage of containers, among all the outbound 
containers, having arrived at a SY until the t* day before the 
arrival of a vessel) x (the length of a period). s$ for t > 0 can 
be evaluated by using {(1 — the cumulative percentage of 
containers, among all the inbound containers, having left the 
terminal until the (t— 1)" day after the departure of a vessel) x 
(the length of a period)}. 

As an example, we illustrate how shy for t < 0 can be 
estimated as follows. In order to evaluate s$ for t < 0, we must 
have the following two data: the operation time for a YC to 
receive an outbound container from an external truck; and the 
percentage of containers, among all the outbound containers, 
arriving at a SY on the t® day before the arrival of the vessel. 
Tab. 4 shows an example of s;.,,, for t< 0. The expected handling 
time for a receiving operation by a YC used in the example is 
1.521 minutes. 


Tab. 4. Example Of Soy fort<0 


t (day) 
Percentage of 


containers [%] 
Shy (minutes) | 0.030 


0.091 | 0.152 


Tab. 5 summarizes the data used for calculating a resource 
profile. Tabs 4 and 6 show the percentage of containers 


arrived at a SY on the days before loading and of containers respectively. Tab. 7 illustrates the resource profile for berthing 
left the terminal on each day after the departure of a vessel, ofa vessel. 
a) consumption 
of resource 


Berth 


Qc 


TR 


YC handling (receiving) 


[] ' 
! ! 
SY i storage (outbound containers) i 
! I 


receivi | loading | period 
| 
' 
1 
L] 
[i 
i 
i 
' 
I 
L] 
[i 


consumption 
b) 
of resource 


Berth berth 


SY storage (inbound iners) 


| unloading | | delivery | period 
Fig. 3. Resource profile for berthing a vessel. a) resource profile for outbound containers, b) resource profile for inbound containers 


Tab. 5. Data used for calculating the resource profile for a vessel in a berth planning 


Length of the vessel plus the allowance between adjacent vessels 
Berthing duration of the vessel 
Number of loading containers for the vessel 
Number of unloading containers for the vessel 
Time for a QC to transfer an outbound container to a slot of the vessel 
Time for a QC to transfer an inbound container to a TR 


Turnaround time for a TR to the travel between the shore and the yard 
Time for a YC to receive an outbound container from an external truck 1.521 min. 
Time for a YC to transfer an inbound container to an external truck 2.242 min. 


Time for a YC to transfer an outbound container to a TR 1.134 min. 
Time for a YC to receive an inbound container from a TR 1.114 min. 
Storage duration of a container during a period (the length of a period) 


Tab. 6. Percentage of inbound containers left the terminal on the n'" day after unloading 


n 1 2 3 4 5 6 
Percentage (%) 34 22 15 14 4 


ute) 


0.076 | 0.091 
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A numerical example is provided to illustrate a resource 
profile for a berth planning. We assume that the length of one 
period is 1 day, and the planning horizon is 14 days. Also, let 
us suppose that a vessel is scheduled to berth at the quay during 
period 7, 01:00~19:00. The total length of the quay is 1,500 m 
and thus the capacity of the berth is 1,500 x 1,440 m-min per 
period. The total number of slots in the SY is 21,000 slots. 
The SY consists of 20 storage blocks in which a storage block 
consists of 5 tiers, 6 rows, and 35 bays. Thus, the capacity of the 
SY becomes 21,000 x 1,440 slot-min per period. The number 
of QCs, TRs, and YCs is 15, 40 and 20 respectively. However, 
a QC has a schedule for preventive maintenances on periods 6, 
7, and 8, and five TRs are expected to enter a maintenance shop 
for such preventive maintenances in periods 12, 13 and 14. The 
capacity of QCs, TRs, and YCs in a period can be evaluated by 
multiplying the total number of equipments with the length of 
a period, namely 21600, 57600 and 28800 machine-minutes, 
respectively. In the periods of 6, 7 and 8, the capacity of QCs 
becomes 20160 machine-minutes. In the periods of 12, 13 and 
14, the capacity of TRs becomes 50400 machine-minutes. 

Tab. 8 shows a part (periods 4 — 10) of an example in 
a capacity plan for the berthing of a vessel in the 7" day. The 
notations ‘CP (capacity)’, ‘RS (reserved)’, ‘AV (available)’, 
and ‘RQ (requirement)’ represent the total capacity of each 
resource during a period, the amount of the resource already 
reserved by other plans, the amount of available resource that 
is CP subtracted by RS, and the amount of the resource that will 
be reserved by the berthing of the vessel under consideration, 
respectively. Values in the ‘RQ’ are calculated by using the 
resource profile noted in Tab. 7 and the number of loading and 
unloading containers for the vessel noted in Tab. 5. 


In Tab. 8, note that the resource requirements for the QCs 
and SY exceed their available capacities in the period of 7. Thus, 
we have to prepare extra QCs and storage spaces for containers 
in the period of 7 or reject the berthing of this vessel. Because 
a QC has a schedule for preventive maintenances in the period 
of 7, changes in a maintenance schedule for the QCs may solve 
these problems. Also, leasing the outside of SY (e.g. off-dock 
container yard) or reducing the dwell time of containers may 
increase in storage spaces for containers. 


SCHEDULING OF QC WORKS (SPLIT) 


Scheduling of QC works determines a schedule for each 
QC to discharge and load containers on a specific deck or in 
a specific hold of a vessel. We assume that a related berth 
plan has already been constructed. Tab. 9 shows the input 
parameters, decision variables, objectives, and constraints for 
scheduling of QC works. 

For constructing a QC work schedule, a stowage plan for 
vessels and yard map in a SY must be provided. A stowage 
plan indicates the location of slots which containers must be 
discharged from or to be loaded in. A yard map shows storage 
locations of containers bound for a vessel. There may be some 
precedence relationships to be satisfied in unloading and loading 
tasks. For example, when a discharging operation is performed 
at a bay in a vessel, operations on a deck must be performed 
before the operations in hold of the same bay start. Also, loading 
operations in a hold must precede the loading operation on 
a deck of the same bay in a vessel. Because QCs travel on the 
same rail, two adjacent QCs must be apart from each other by 
at least a certain distance for them to simultaneously perform 


Tab. 8. Example of a capacity plan in a berth planning 
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21,600 


43,185 
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36,075 


14,415 
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21,360 
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28,800 


28,800 


28,800 


28,800 


28,800 


28,800 


28,800 


24,900 


25,380 


24,960 


23,610 


22,140 


22,020 


21,540 


3,900 


3,420 


3,840 


5,190 


6,660 


6,780 


7,260 


82 


99 


534 


1,236 


427 


276 


188 


30,240,000 


30,240,000 


30,240,000 


30,240,000 


30,240,000 


30,240,000 


30,240,000 


28,800,000 


29,160,000 


29,088,000 


28,800,000 


28,080,000 


27,720,000 


27,360,000 


1,440,000 


1,080,000 


1,152,000 


1,440,000 


2,160,000 


2,520,000 


2,880,000 


56 


178,848 
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272,160 


777,600 


1,584,000 


806,400 


532,224 


354,816 


their operations without any interference. Planners attempt 
to construct in a way of minimizing the makespan in ship 
operations. 


Tab. 9. Definition of the scheduling problem of QC works 


Stowage plan of a vessel 


Input Available time window of each QC 


parameters | Yard map 
Resource profiles 


Decision 


E Work schedule for QCs assigned to a vessel 


Objectives | To minimize a make-span in ship operations 


Precedence relationships in operations 
Constraints | Interference among QCs 
Availability of resources 


Like in the berth planning, resource requirements depend 
on the number of unloading and loading containers and are 
evaluated by using a load profile for the scheduling of QC 
works. Requirements for each resource can be represented as 
shown in Figs 4-(a) and (b) that illustrate the resource profiles 
applied in inbound and outbound flows, respectively. 

a) consumption 
of resource 


period 


consumption 
of resource 


b) 


YC 


TA 


1 |] 
loading „| period 
Fig. 4. Resource profile for unloading and loading operations in a vessel. 


a) resource profile for inbound containers, 
b) resource profile for outbound containers 


During the process of the scheduling of QCs in advanced 
container terminals, blocks where inbound containers are 
unloaded and outbound containers are picked from each ship- 
bay are simultaneously determined. The availability of the 
resources in each corresponding block must be checked before 
the decision on the blocks is fixed. The resources in each block 
are Y Cs, TAs, and SY of each block. However, because the fleet 
of TRs may be operated in a pool for all QCs, the capacity of 
all TRs in a terminal are compared with the required amount 
of TRs during the ship operation. 


Let sĝy, and sĝy, be the unit amount of the YC capacity 
required for transferring a container for unloading and loading 
operations that is required at the period when the ship operation 
occurs, respectively. so ,, and Soa, can be defined in the same 
way as previously mentioned ways. A resource profile can be 
evaluated as follows: s{. can be evaluated by using {(the time 
for a QC to transfer an outbound container to a slot of a vessel) 
x (the percentage of loading containers among all containers 
for a QC) + (the time for a QC to transfer an inbound container 
to a TR) x (the percentage of unloading containers among all 
containers for a QC)}. Because TRs are required in loading and 
the unloading operations, s2, can be evaluated in the similar 
way as Sic. Soy, and sĝy, can be evaluated by using the time for 
a YC to receive an inbound container from a TR and time for 
a YC to transfer an outbound container to a TR, respectively. 
Because TAs are consumed at the same period when YCs are 
required, s¢,, and sĝą can be evaluated in the similar way as 
Soy, and sy, respectively. 

~ Tab. 10 summarizes data used for calculating the resource 
profile. For the transfer time of a QC, cycle time of a TR, 
and handling time of a YC, we used those in Tab. 5. Tab. 11 
illustrates a resource profile in ship operations by QCs. 


Tab. 10. Data used for calculating a resource profile 
in the scheduling of QC works 


Number of loading containers for a QC | 15 during period 3 


Number of unloading containers 
for a QC 
Time for a TR to transfer a loading 
container ina TA 
Time for a TR to transfer an unloading 
container ina TA 


16 during period 3 


3 min 


3 min 


Tab. 11. Resource profile in ship operations by OCs (So ) (unit: minute) 


KA eee eee 


aS ee 


We will introduce a numerical example to illustrate the 
resource profile applied to the scheduling of QC works. We 
assume that the length of one period is | hour, and a planning 
horizon is 12 hours. Also, let us suppose that a QC is scheduled 
to handle the containers during the period of 3. The capacity 
of a QC is 60 machine-minutes, the capacity of TRs is 2,400 
machine-minutes, and the capacity of YCs in a storage block 
is 60 machine-minutes (refer to Section 3.1). Assuming that 
a storage block has 9 TAs and average occupation time ofa TR 
in the TA is 3 minutes, the capacity of TAs in a storage block 
for a period can be obtained by multiplying the number of TAs 
in a storage block with the maximum number of TRs passing 
the TA for a period, namely, 180 machine-minutes. 

The containers to be discharged or loaded by the QC will be 
located or are located in four different storage blocks as shown 
in Tab. 12. Tab. 13 shows an example of a capacity plan for 
the scheduling of QC works. Entries of ‘RQ (requirement)’ of 
each resource are calculated by using the resource profile in 
Tab. 11, number of loading and unloading containers for the 
QC in Tab. 10, and storage distribution of the containers for 
the QC in Tab. 12. 


Tab. 12. Storage distribution of the containers for the OC 


Storage block 


Number of unloading containers 
Number of loading containers 
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Tab. 13. Example of a capacity plan for the scheduling of QC works 


In Tab. 13 note that the resource requirements for the YCs 
in the storage block of 1 exceed the available capacity in the 
period of 3. Therefore, we have to deploy additional YCs in the 
period of 3. If other storage blocks can be used for unloading 
and loading containers in the period of 3, then the planner may 
modify the original schedule. Of course, the availability of the 
YCs and TAs must be checked again after the modification. 
Unfortunately, if a terminal cannot provide additional YCs in 
the period of 3, a certain part of work schedule in the period 
of period 3 must be moved to other periods. 


YARD PLANNING 


There are different types of yard plans depending on the 
handling facilities in a yard. However, most of yard plans specify, 
at least, the number of inbound containers to be discharged from 
each vessel and then stacked at each block, and the number of 
outbound containers bound for each vessel to be stacked at each 
block. Some yard plans may provide more detailed decisions 
on the storage, such as plans on storage layouts of outbound 
containers in each block. This paper assumes that a yard plan 
is constructed before containers bound for or discharged form 
a vessel start to arrive at a specific yard. 

The decision on the exact storage location of an individual 
container is made in real time whenever a container arrives at 
a yard. Tab. 14 shows input parameters, decision variables, 
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objectives, and constraints for a yard planning. For making 
a good yard plan, planners attempt to balance the workload 
in storage blocks and to minimize the travel distance for TRs 
to transport the containers between their storage blocks and 
berthing locations of the corresponding vessels. 


Tab. 14. Definition of the problem in a yard planning 


Number of outbound containers of each group 
that are bound for a vessel and will arrive at 
a terminal during each period 

Number of inbound containers discharged 
from a vessel during each period 

Resource profiles 

Storage amounts for each group of containers 
to be stored in each storage block 

To balance workloads among storage blocks 
To minimize travel distances for TRs to 
transport containers between storage blocks 
and vessel berthing locations 

Availability of resources 

Storage capacity of each storage block 
Handling capacity of each storage block 


Input 
parameters 


Decision 
variables 


Objectives 


Constraints 


The requirement of each resource can be determined 
as a time-phased way with respect to the arrival time of an 


a) consumption 
of resource 


1 
1 
TR ! 
1 
YC handling (receiving) 
[j 
TA transfer (receiving) 
[j 
1 
SY i 
[] 
b) consumption 
of resource ; ; 
= [j 
ri po ' 
™ (berth to SY) 
[j 
YC (unloading) 
' 
TA | (unloading) | 
SY i 
[] I t 
unloadii 


Fig. 5. Resource profile for inbound and outbound containers. a) resource profile for outbound containers, b) resource profile for inbound containers 


outbound container or the discharging time of an inbound 
container. Figs 5-(a) and (b) illustrate a resource profile for 
outbound and inbound flows, respectively. As in the case of 
the scheduling of QC works, it is necessary to consider the 
YCs, TAs and storage spaces by each storage block and all 
TRs in a terminal. Decision variables in a yard planning are 
the number of outbound or inbound containers to be stacked 
at each block in each time period. 

A resource profile for outbound containers can be evaluated 
as follows: Sbg for t > 0, can be evaluated by using the 
turnaround time for a TR to the travel between the shore and 
the yard. spy for t> 0 can be evaluated by (the transfer time for 
a YC to receive an outbound container) x (the percentage of 
containers, among all the outbound containers, arriving at the 
SY on the (t + 1)" period from the starting period of arrivals). 
spy for t> 0 can be evaluated by using the time for a YC to 
transfer an outbound container from a TR. Because TAs are 
required at the same period when YCs are required, Spa for 
t > 0 can be evaluated in the similar way as spy for t > 0. Sbs 
for t > 0 can be evaluated by (the cumulative percentage of 
containers, among all the outbound containers, arriving at the 
SY until (t + 1)" period from the starting period of arrivals) x 
(the length of a period). 

Second, a resource profile for inbound containers can be 
explained as follows: Spr for t = 0 can be evaluated by using 
the time for a TR to the travel between shore and yard. spy for 
t = 0 can be evaluated by using the time for a YC to receive an 
inbound container from a TR. spy for t > 0 can be evaluated by 
using (the time for a YC to transfer an inbound container to an 
external truck) x (the percentage of containers, among all the 
inbound containers, leaving the terminal on the t" period after 
unloading). Because TAs are consumed at the same period 
when YCs are required, sba for t > 0 can be evaluated in the 
similar way as Spy for t > 0. Shs for t = 0 can be evaluated by 
using the storage duration of the container for a period. Sbs, for 
t > 0, can be evaluated by {(1 — the cumulative percentage of 


containers, among all the inbound containers, having left the 
terminal until the (t— 1)" period after unloading) x (the length 
of a period)}. 

Tab. 15 summarizes data used for calculating the resource 
profile. For the cycle time of a TR, handling time of a YC, 
and transfer time of a TR, we used ones in Tabs 6 and 11. 
Tabs 16 and 17 illustrate the resource profile for receiving and 
discharging of containers, respectively. 


Tab. 15. Data used for calculating a resource profile in a yard planning 


Number of receiving containers 
for a vessel 


Number of unloading containers 


for a vessel >69 


Time for an external truck to transfer 


ie ie 3 minutes 
a receiving container in a TA 


Time for an external truck to transfer 


i Sep 3 minutes 
a delivery container in a TA 


24 hours or 
1,440 minutes 


Storage duration of the container 
for a period 


Tab. 16. Resource profile for receiving 
an outbound container (s:,,) (unit: minute) 
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Tab. 17. Resource profile for discharging 
an inbound container (51) (unit: minute) 


The following introduces a numerical example to illustrate 
a resource profile in a yard planning. We assume that the length 
of one period is | day and the planning horizon is 9 days. Let us 
suppose that a yard is scheduled to receive outbound containers 
from period 2 to 7 and to discharge inbound containers during 
the period of 2. The capacity of a storage block in the SY 
becomes 1,050 x 1,440 slot-minutes — a storage block with 
5 tiers, 6 rows, and 35 bays. The capacity of TRs is 57,600 
machine-minutes, the capacity of YCs in a storage block is 
1,440 machine-minutes, and the capacity of TAs in a storage 
block is 4,320 machine-minutes. 

Planners decide that inbound containers are stacked in 3 
storage blocks, and outbound containers are stacked in 5 storage 
blocks as shown in Tab. 18. Tab. 19 illustrates an example of 
a capacity plan in a yard planning. Entries of ‘RQ (requirement)’ 
of each resource correspond to the amounts of the resources that 
are required for the storage in Tab. 18. The values in ‘RQ’ are 
calculated by using the resource profile in Tabs 16 and 17, the 
number of inbound and outbound containers that are planned to 
be stacked at the block 1 and 3 as shown in Tab. 18. 


Tab. 18. Storage plans for inbound and outbound containers 


Storage 1 
block 


Number 
of inbound | 220 


containers 


Number 
of outbound 
containers 


In Tab. 19 note that the resource requirements for the YCs 
and TAs in the storage block 1 exceed their available capacities 
in the period of 2. Therefore, we have to prepare extra YCs and 
TAs in other periods. 


CONCLUSIONS 


— For constructing an efficient operational plan in container 
terminals, a large number of factors must be considered for 
decision-making. Although various planning activities are 
mutually related with each other, they have been treated as 
independent decision-making processes. Furthermore, they 
often share the same resource that they have to compete 
with each other to secure. 

— This study proposed a unified framework for integrating 
various planning activities in container terminals and 
defined decision-making problems in each planning 


Tab. 19. Example of a capacity plan in a yard planning 


a) TR 


57,600 57,600 57,600 57,600 


57,600 57,600 57,600 57,600 57,600 


47,160 43,200 36,540 36,240 


36,075 36,015 32,610 32,610 


10,440 14,400 21,060 21,360 


21,525 21,585 24,990 24,990 57,600 


5,600 = 


4,320 


3,700 


620 


224 


92 73 


1,512,000 


1,512,000 


1,512,000 | 1,512,000 


1,512,000 


1,512,000 | 1,512,000 | 1,512,000 


1,157,000 


1,140,000 


1,104,000 | 1,103,000 


1,102,000 


1,101,000 | 1,100,000 | 1,099,000 


355,000 


372,000 


408,000} 409,000 


410,000} 411,000} 412,000) 413,000 


—| 316,800 


316,800] 209,088 
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139,392 91,872 47,520 12,672 


Tab. 19. Example of a capacity plan in a yard planning 


c) YC, TA, and SY in the storage block 3 


4,320 


3,900 


420 


161 


1,512,000 | 1,512,000 | 1,512,000 | 1,512,000 


1,512,000 


1,512,000 | 1,512,000 


1,512,000 


1,512,000 


1,139,000 | 1,140,280 | 1,102,000 | 1,101,000 


1,100,000 


1,090,000 | 1,080,000 


1,070,000 = 


373,000 371,720] 410,000} 411,000 


412,000 


422,000} 432,000} 442,000) 1,512,000 


— 2,592 9,072 16,848 


activity. Input parameters, decision variables, objectives, 
constraints, time buckets, and planning horizon for each 
decision activity were identified. The concept of a resource 
profile was suggested that should be utilized to check 
the feasibility of a plan with respect to the constraint on 
the resource availability. Numerical examples were also 
provided to illustrate capacity planning procedures for 
a berth planning, scheduling of quay crane operations, and 
yard planning. 

— The concept proposed in this study may be utilized for 
developing operational software for container terminals. 
Also, methods for optimizing various operational decisions 
by utilizing the concepts proposed in this study may be 
developed in future studies. 
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Occurrence probability of maximum sea levels 
in Polish ports of Baltic Sea coast 


Bernard Wisniewski, Prof. 
Tomasz Wolski, Ph. D. 
University of Szczecin 


ABSTRACT 


In this work long-term probability of occurrence of maximum sea levels in some points of Polish Baltic 
Sea coast, was determined. Use was made of multi-year series of measurement data on maximum yearly 
sea levels, and their probability distributions were determined. To the analysis Gumbel s distribution and 
Pearson distribution of 3rd type as well as quantile methods and the highest credibility method, were applied. 
Kolmogorov test was used to examine conformity of the theoretical distributions with real random variable 
distribution. As results from the analysis, the highest sea levels of 1000- year return period can be expected 
in Polish ports of the west part of the coast , i.e. Kolobrzeg (750, 2 cm, i.e. 2,5 m above the average sea 
level) and Swinoujscie (723,6 cm) . Lower sea levels of the same return period can be expected in Ustka 
(720,2 cm), Wladyslawowo (709,7 cm) and Gdansk (716, 7 cm), respectively. 


Keywords: sea levels; probability; Polish Baltic Sea coast 


INTRODUCTION 


This work is aimed at determination of long-term 
probabilistic prediction of maximum yearly sea levels in five 
Polish Baltic Sea coast observation stations: Swinoujscie, 
Kolobrzeg, Ustka, Wladyslawowo and Gdansk. Knowledge 
of sea levels of a given occurrence probability is necessary in 
designing the marine hydro-engineering buildings as well as 
in determining the characteristics of storm swellings. 

The most important is to know high water level stages which 
can happen once a determined number of years, e.g. once 50, 
100, 150 or 200 years. This period is called the return period T. 
It is expressed as follows: 


100 
= — [years] (1) 
p% 
where: 
T — mean return period 
p — probability [%.]. 


Empirical probability is defined by Weibull ’s function [1]: 


p(m,N)% = 100m (2) 
N+l 
or Czegodajew formula used in hydro-engineering [7]: 
-025 
nim. = n G) 
N+0.50 
where: 
m — successive term of distribution series 
N — size of series. 
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To estimate a probable water level, available measurement 
data series from previous years should be used. Credibility 
of such estimation depends both on length of water level 
observation series and representativeness of a chosen 
theoretical probability distribution function. In determining 
probability of the yearly highest sea levels a random sample 
covering a few dozen years of measurements is usually at 
one’s disposal. Hence for Kolobrzeg the observed water 
levels from the years of 1867 + 2006 were at one’s disposal, 
for Swinoujscie — from the years of 1901 + 2006, for Gdansk 
— from the years of 1886 + 2006, for Ustka — from the years 
of 1946 + 2006 and for Wladyslawowo — from the years of 
1947 + 2006. 

Calculations of exceedance probability of water levels 
consist in appropriate selection of theoretical probability 
distributions and next in assumption of estimation methods of 
parameters of a selected distribution, with the use of statistical 
data. Confirmation of the fact whether the distribution has been 
properly selected can be obtained by using the tests of goodness 
of fit, proposed in the literature. Calculation results are given 
together with a confidence interval within which real variable 
value can be found with a given probability. 

In probability theory the random variable X is such quantity 
which - due to random coincidence of various factors — can 
achieve different numerical values with a given probability. 
In this work the series of yearly maximum water levels 
observed at the stations of Swinoujscie, Kolobrzeg, Ustka, 
Wladyslawowo and Gdansk, are taken as the random variable X. 
A function which determines how large is occurrence 


probability of the event that the random variable X takes one 
of the values contained within the numerical interval S on the 
x - variable axis is called the probability distribution of the 
random variable X. 
In hydrology and oceanology the distributions: 
— log-normal distribution (Gauss — Laplace) 
— Pearson’s distribution of type II 
— Debski’s distribution 
— Gumbel’s distribution of extreme values (Fisher-Tippett, 
type I) are most often used [1, 6, 8]. 


Every theoretical distribution density function has a number 
of parameters g,,...,g,, on which its graphical form depends. 
One of the crucial tasks is, on the basis of random sample, 
to estimate values of the paramters g.,...,g,, in such a way 
as random features of general population to be represented 
as best as possible by the function f(x, g,,...,g,). In statistics 
to this end serve the methods of estimation of probability 
distribution parameters. Usefulness of particular estimation 
methods is dependent first of all on random sample size as 
a well as on a required degree of approximation of estimated 
values to theoretical ones. To large random samples simpler 
methods which provide lower estimation accuracy, can be 
applied [7]. 

In practice it is not possible to fully confirm how far 
a theoretical distribution fits actual distribution of random 
variable, i.e. that empirically elaborated on the basis of 
statistical observation series. In such situation consistency 
of theoretical distribution and observation results should be 
examined. To this end serve statistical consistency tests which 
determine occurrence probability of differences between 
theoretical distribution and empirical one. The A point test of 
Kolmogorow and the y’ linear test of Pearson belong to the 
most often used [1, 8]. 

A few calculation methods of water levels with a given 
exceedance probability can be distinguished depending on 
an assumed type of probability distribution curve as well as 
estimation method of parameters. Among the most frequently 
applied the following methods can be numbered: 

— quantile method (of position characteristics) 
— maximum likelihood method and 
— method of statitstical moments. 


In this work Pearson’s distribution and Gumbel’s one was 
used and from among the applied methods the quantile one and 
that of highest likelihood was selected as recommended for the 
determining of occurrence probability of extreme hydrological 
phenomena such as water flow rates and water levels [6], [7]. 


Quantile method 


The method is based on determined specific values of the 
variable X called quantiles. The p — order quantile is called such 
value of the variable x for which the exceedance probability 
is equal to p. Distribution quantiles are compared to sample 
quantiles. To determine quantiles from random sample its 
elements should be ordered according to their values to form 
the so called distribution series. The following quantiles are 
usually determined: 
X <,, — a value corresponding to the middle of series, 
equivalent of median 
= X y X 100p WO symmetrical quantiles, equally distant from 
the middle of series, where: p = 10% or p = 5% is selected 
most frequently 

— X jo, X oy — two quantiles equivalent to lower and upper 
limit values of random variable. 


If distribution series is long nad regular then values of 
symmetrical and middle quantiles can be determined by using 
interpolation between terms of the series. In the case of short 
and irregular series the searched - for quantiles are read from 
the equalized diagram of summed-up frequencies. However 
the values x oo, aNd X „are usually determined subjectively 
as they are located beyond the range of values of sample 
elements [1, 8]. 


Maximum likelihood method 


The basic notion of the method is the likelihood function: 


L=£ (6), 8)5--8p)* IE Spor -Bq) *- E Spoon) = 


x (4) 
=] PO j.81.--8,) 
1 


where: 

f(x, g,,...g,) — probability density function of the random 
variable x, in a given mathematical form and 
having unknown distribution parameters 


£,(2,.8,--8,) 


The function can be approximately interpreted as probability 
of obtaining the random sample Z (x) exactly the same as that 
really observed. For calculation reasons to use the logarithm 
of likelihood function is more convenient: 


InL =Inf(x,,g),....8,)+..-Inf(«y, 81)... 82) = 


N 
= X in FR gogn 


Estimation of distribution parameters consits in finding 
such system of the values g,, g,,..-g, at which probability of 
observing the given random sample Z (xi) is the greatest. The 
discussed method is equivalent to searching for maximum 
value of the function L respective to g,value, i.e. to solving 
the set of equations: 


(5) 


=0 (6) 


The maximum likelihood method involves great troubles 
in the case when one of the estimated parameters is lower or 
upper variability range limitation of the random variable X. 
In such cases to use other methods is proposed for estimating 
the limitations [1, 8]. 


APPLIED METHODS 


Determination of theoretical probable maximum 
sea levels by means of the quantile method with 
the use of Gumbel’s distribution 


The Gumbel’s distribution density function is based on 
statistical distributions of extreme values which occur in certain 
larger sets of values. For instance it can be maximum values of 
sea levels considered in this paper. The Gumbel’s distribution 
density function is double exponential and described by the 
formula [21: 


A 


x—b 


x—b 
—-e 


toate- (7) 
a 
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where: 

A — scale parameter (it determines dispersion of the 

7 distribution along x-axis) 

k — location parameter (it determines location of the 
distribution on x-axis) 

e — Napierian base. 


Essence of estimation of the assumed distribution against 
given measurement data is to determine estimators of the 
distribution parameters â and F. 

After finding logarithms and simplifying the above given 
formula the following was obtained: 


Xo, =F -1/ a * Inf-In(1-fx))] (8) 


To estimate the parameters 4 and F, formulas for the three 
quantiles: upper one X, middle one X, and lower one X, 
were preliminarily determined. Here the following values were 
assumed: p = 5% and 1-p = 95%. 

In compliance with the formula (8) for X „ the following 
was achieved [7]: 


X „= Ê -1/a(-2.9702) (9) 
X soy, =Ê -1/A(-0.36651) (10) 
X yo, = Ê -1/a(1.097189) (11) 


Next the variability measure was determined form the 
formula: 


V=(X,-X,_y/2 (12) 


By substituting relevant expressions from Eqs (9) and 
(11) for X, and X, „ respectively, the following was obtained: 


V = 4.0672/24 — A = 2.0336/V (13) 


Then, on the basis of distribution series of maximum water 
levels for each of the considered stations the water levels were 
read or interpolated by using the empirical probability from the 
formula (2), for X, and X,_, at p = 5%. The obtained values 
were used for determining the distribution parameters â and F 
from Eqs. (12), (13) and (10), respectively. The so achieved 
distribution parameters â and F were introduced to Eq. (8) 
for X „ and the theoretical sea levels for selected probability 
quantiles were calculated for each of the observation stations 
(Tab. 1, 2,3, 4 and 5) 


Tab. 1. Theoretical maximum sea levels and probability of their occurrence 
at Swinoujscie observation station in the period of 1901+2006 years 
Sea level [cm] 


Gumbel’s Person’s 
distribution distribution 


vo 
= 


Q 
= 


T (years) 


method 


OF 


OF 


likelihood 
likelihood 


Maximum 


Sea level [cm] 


Gumbel’s Person’s 
distribution distribution 


~ 
A 
= 
Ss 
vo 

2 

= 


Maximum 
likelihood 
method 


Tab. 2. Theoretical maximum sea levels and probability of their occurrence 
at Kolobrzeg observation station in the period of 1867+2006 years 
Sea level [cm] 


Gumbel’s 
distribution 


Person’s 
distribution 


T (years) 


Maximum 
likelihood 
method 


75% 578.5 576.4 


90% 
95% 


567.2 
561.3 


564.2 
555.0 


99% 
100% 


551.7 
537.6 


546.1 
537.6 


Tab. 3. Theoretical maximum sea levels and probability of their occurrence 
at Ustka observation station in the period of 1946=2006 years 


Sea level [cm] 


Gumbel’s 
distribution 


Person’s 
distribution 


T (years) 


likelihood 
likelihood 
method 


Maximum 
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Tab. 4. Theoretical maximum sea levels and probability of their occurrence 
at Wladyslawowo observation station in the period of 1947+2006 years 


Sea level [cm] 


Gumbel’s 
distribution 


Person’s 
distribution 


T (years) 


method 
likelihood 
method 


Maximum 


Maximum 
likelihood 


Vx! 


(14) 


where: 

a, g, 4 — distribution parameters which should satisfy the 
following conditions: x > e (lower limitation of the 
distribution), a > 0, à > 0 

Tà) — gamma function of the variable À. 

In this work the simplified form of the function which 
determines Pearson’s 3rd type probability distribution, was 
used [5, 6]: 

X *[1 + PẸ, s)* C] 


p% X so% (15) 

where: 

X œ — value of the middle quantile, 
— variability coefficient, 

o (p.s) — function of exceedance probability and skewness 


(asymmetry). 


The following are parameters of the distribution: 


Tab. 5. Theoretical maximum sea levels and probability of their occurrence 
at Gdansk observation station in the period of 1886+2006 years 
Sea level [cm] 


Gumbel’s 
distribution 


Person’s 
distribution 


T (years) 


Maximum 
likelihood 
method 


Determination of theoretical probable maximum 
sea levels by means of the quantile method with 
the use of Pearson’s distribution of type IIT 


In hydrology Pearson’s distribution is used most often. Its 
type III has found the widest application to equalization of 
series of numerical features of hydrological phenomena. The 
density function of Pearson’s 3rd type distribution is of the 
following form [5, 6]: 


the variability measure 
v= Xion - Xooy,)/ 2 
— the variability coefficient 
C= VK ox (16) 


— the skewness coefficient s determined from an auxiliary 
table depending on value of the expression: 


C, = (C, * Xo); 50% Xoo) (17) 


Like in the Gumbel’s distribution, on the basis of the 
distribution series of maximum water levels for each of the 
observation stations the water levels were read or interpolated by 
using the empirical probability from Weibull’s formula (2), for 
X and X, » at p= 10% (the quantiles: X 0 X 50% X o0 X100% 
The obtained values of quantiles were used for determining the 
distribution parameters: V, C, and s from Eqs. (12) and (16) 
as well as (17) and the auxiliary table, respectively [5]. The 
calculated distribution parameters and the function ®(p,s) read 
from the above mentioned auxiliary table, were introduced to the 
successive formula, (15), for X , and the theoretical sea levels 
for the selected probability quantiles for each of the observation 
stations were calculated (Tabs 1, 2, 3, 4 and 5). 


Determination of theoretical probable maximum 
sea levels by means of the maximum likelihood 
method with the use of Pearson’s distribution 


of type IIT 


In this method the form of Pearson’s distribution function, 
proposed by Foster, was applied [6, 1]: 


X,=et 1/a*tp (18) 
where: 

e — lower limitation of the distribution, equivalent of the 
quantile X o0% 

value ofexceedance probability function and parameter A, 


read from the auxiliary tables [6, 9] 


t= 
P 


a — the second parameter of the distribution, determined from 
the formula: 
a= A/x,. (19) 
where: 


x „ — arithmetic mean of the values x,, where x, - the substitute 
variable determined from the formula: 
x; = x max- £ fori =1, 2,3,...N (20) 
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where: 
x, — maximum sea levels in distribution series as well as the 
parameter A determined from the formula: 


R= 1/(4*A) * [1 + (1 + 4/3 *A)0.5] QD) 
where: 
A — function defined by the formula: 
A= In(x,) - (nx), (22) 
where: 
lnx, — logarithm taken from the arithmetic mean of the 


values x, 
(lnx), — arithmetic mean of the values Inx. 


In the maximum likelihood method, on the basis of the 
distribution series of maximum water levels, for each of the 
observation stations the substitute variable x,, acc. Eq. (20), and 
its arithmetic mean x „ as well as the logarithm of the mean, 
lnx „and the arithmetic mean of the logarithm values, (Inx)sr, 
then also the functions A, acc. Eq. (22), and the parameter ë, 
acc. Eq. (21), which was used to calculate the second parameter 
of Pearson’s distribution, a, acc. Eq. (19), were determined. 
The lower limitation of the distribution, l, was determined 
from the equation of the trend-line of maximum sea levels, 
depicted on the quantile - quantile diagram (by using Statistica 
software). 

The functions t were read from the auxiliary tables [6, 9]. 
Having all the Pearson’s distribution parameters one made use 
of the Foster’s formula, acc. Eq. (18), and calculated theoretical 
sea levels for selected probability quantiles for each of the 
observation stations (Tab. 1, 2, 3, 4 and 5). 


Test of consistency of the assumed theoretical 
probability distribution with empirical 
distribution 


In this work consistency of the assumed theoretical 
distribution with empirical one (observation series of sea 
levels) was examined with the use of Kolmogorow consistency 
test. The testing consists in checking the following condition 
[1, 8]: 

Dmax LP/m, N/% — p%]< A, WN (23) 
where: 
p/m, N/ % — empirical occurrence probability of m-th term of 
distributive series 


p% — theoretical occurrence probability of sea 
level value corresponding with m-th term of 
distributive series 

he — critical value of the Kolmogorow distribution 
(at a = 5% — kh, = 136, at a = 1% — d, = 163) 

N — size of distributive series. 


The condition of the Kolmogorow test was checked for all 
the methods of determination of theoretical probable maximum 
sea levels for each of the observation stations. Results of the 
test are given in Tab. 6. 


Calculation of estimation error 
and confidence interval 


Confidence interval is a range within which real value of 
quantile will occur with the assumed probability P (confidence 
level). In this work the confidence interval limit x was 
determined on the level P, = 68.3% as well as P, = 95% by 
means of the following formula [1, 4]: 

for the upper confidence interval: 


Xp =Xpt ta* © (Xp) (24) 
for the lower confidence interval: 
Xp = Xp ta* 0 (Xp) (25) 


for the confidence level = 68.3% 

for the confidence level = 95% 

o — estimation error determined according to the 
following formula: 


Ox F(s,p)*C,* xsov AN (26) 


where: 
F(s,p) — function of probability and assymetry measure, read 
from the auxiliary tables [4] 


The determined confidence interval limits for the confidence 
levels P, = 68.3% and P, = 95% are given in Tab. 7. 


ESTIMATION RESULTS OF THEORETICAL 
MAXIMUM WATER LEVELS FOR 
PARTICULAR OBSERVATION STATIONS 
AND OCCURRENCE PROBABILITY OF 
THE LEVELS 


In this work the sea levels with assumed occurrence 
probability were determined on the basis of Gumbel’s 
distribution and Pearson’s distribution by applying both 
the quantile and maximum likelihood methods.To check 
consistency of the assumed theoretical distributions with 
empirical ones the Kolmogorow test was used (Tab. 6). Results 
of the test calculations do not lead to rejection of the in-advance- 
assumed hypothesis on consistency of the distributions. 

The best adjusted to long-term observation series are the 
theoretical maximum sea levels and their occurence probability 
determined by means of the Gumbel’s distribution and maximum 
likelihood method, that is confirmed by description given in 
the chapter „Applied methods” and the diagrams presented 
in Fig. 1 through 5. The so-calculated highest theoretical sea 


Tab. 6. Results of Kolmogorow test for various methods of determination of theoretical sea levels at particular observation stations 


Ap NN 


Observation stations oroS i 


Swinoujscie (1901+2006) 


A NN 
for a = 5% 


Gumbel’s 
distr., 
method of 
maximum 
likelihood 


Pearson‘s 
distr., 
method of 
maximum 
likelihood 


Gumbel’s 
distr., 
quantile 
method 


Pearson’s 
distr., 
quantile 
method 


Kolobrzeg (1867+2006) 


Ustka (1946+2006) 


Wladyslawowo (1947+2006) 


Gdansk (1886+2006) 
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Tab. 7. Limits of confidence intervals, [cm], of maximum yearly sea levels at Swinoujscie, Kolobrzeg, Ustka, Wladyslawowo and Gdansk 
for the confidence levels: Pa = 68.3% and Pa = 95.0% 


P(%) 


Swinoujscie 


Kolobrzeg 


Ustka 


Wladyslawowo 


Gdansk 


levels which can occur once a 1000 years, 200 years, 100 years are located in river estuaries and water level indicators installed 
and 50 years are given in detail in Tab. 1 through 5. From the in them are about 1 km distant from sea coast line. For instance 
results it can be concluded that the highest of them will occur in the port of Kolobrzeg the water level of 750.2 cm, i.e. 2.5 m 
and did really occur in Polish west - coast ports, i.e. Kolobrzeg above the average mean sea level, should be expected once 
and Swinoujscie. The lowest of them will occur in the port of a 1000 years (zero points of Polish water level indicators are 
Wladyslawowo which belongs to the Polish middle coast of based on the Amsterdam’s zero-level equal to -500 N.N.). 
Baltic Sea and is directly exposed to sea. The remaining ports Swinoujscie (with 723.6 cm), Ustka (with 720.2 cm), 

Swinoujscie, Gumbel Distribution, Maximum Likelihood Method Ustka, Gumbel Distribution, Maximum Likelihood Method 
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Fig. 1. Occurrence probability of maximum yearly sea levels Fig. 3. Occurrence probability of maximum yearly sea levels 
at Swinoujscie in the period of 1901+2006 years (Gumbel 5 distribution, at Ustka in the period of 1946+2006 years (Gumbel 5 distribution, 
maximum likelihood method) maximum likelihood method) 
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Fig. 2. Occurrence probability of maximum yearly sea levels Fig. 4. Occurrence probability of maximum yearly sea levels 
at Kolobrzeg in the period of 1867+2006 years (Gumbel 5 distribution, at Wladyslawowo in the period of 19472006 years (Gumbel s distribution, 
maximum likelihood method) maximum likelihood method) 
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Gdansk (with 716.7 cm) and Wladyslawowo (with 709.7 cm) 
are the successive ports in which the maximum thousand-year 
levels can be expected. Other occurrence probabilities can 
be also considered. Graphical representation of the relations 
is presented in Fig. | through 5 for particular Polish ports of 
Baltic Sea coast. 

Gdansk, Gumbel Distribution, Maximum Likelihood Method 
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Fig. 5. Occurrence probability of maximum yearly sea levels 
at Gdansk in the period of 1886+2006 years (Gumbel 8 distribution, 
maximum likelihood method) 


DISCUSSION OF THE RESULTS 


As already mentioned, due to its characteristics Gumbel’s 
distribution has been considered the closest to empirical one. 
Its advantage results from lack of lower limitation of variability 
range of the random variable X, among estimated parameters. 
Estimation of lower or upper limitation of Pearson’s distribution 
(X 100% OF X gy, faces serious difficulty due to a subjective way 
of its determination [1], [8]. 

On the basis of comparison of both methods of estimation 
and determination of distribution parameters it can be concluded 
that the maximum likelihood method makes it possible to 
obtain more precise results as in this method - in contrast 
to the quantile method — to use interpolation of distribution 
parameters is not necessary. 

In Polish subject-matter literature the methods of 
determination of probability of extreme sea levels were 
described in the publications of Wroblewski, A. [10, 11, 
12, 13, 14] and Massel, S. [7]. Prediction of Baltic Sea 
extreme levels was also considered by Jednoral, T., Sztobryn, 
M. and Milkowska, M. [3]. This work presents the methods 
in question in the most complete way as it takes into account 
calculation approaches used by different authors. In this work, 
the characteristics of theoretical maximum sea levels and their 
occurrence probability in five reperesentative Polish ports of 
Baltic Sea coast, were considered. The obtained results are 
based on the longest observation series. From the point of view 
of the applied methods the most important have been so far two 
publications by Wroblewski [11, 13] with which to compare 
results of this work is deemed purposeful. In Tab. 8 can be 
found the comparison of results of the theoretical sea levels, 
determined by Wroblewski for selected quantiles, with those 
obtained by these authors for Gdansk observation station. 

From the comparison of the results of the theoretical 
sea levels for Gdansk observation station, obtained by 
these authors by using the maximum likelihood method and 
Gumbel’s and Pearson’s distributions, with those determined 
by Wroblewski [11, 13], the difference from several up to 
a dozen or so centimeters for different quantiles, can be 
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observed. The differences are mainly caused by a different way 
of determination of the lower limitation of the distribution, 1. 
Wroblewski, determining the lower limitation of the distribution, 
made use of the maximum likelihood method as well as the 
tables of modal values. In this work the equation of the trend- 
line of maximum sea levels, depicted on the quantile-quantile 
diagram obtained from Statistica software, was applied. The 
other cause of the differences in determining the theoretical sea 
levels is length of observation series. In this work the series 
for Swinoujscie, Kolobrzeg and Gdansk are by 37 years longer 
than that used by Wroblewski [11, 13] 


Tab. 8. Theoretical maximum sea levels and probability of their occurrence 
at Gdansk observation station, determined by means of both Gumbel s and 
Pearson 5 distributions with the use of the maximum likelihood method, 
according to Wroblewski [11, 13] and these authors, respectively 


Sea levels [cm] 


T (years) 
Gumbel distr. 
Pearson distr. 
Gumbel distr. 
Pearson distr. 


Results 
by these authors 


ults 
blewski 


CONCLUSIONS 


— In this work sea levels for a given occurrence probability, 
based both on Gumbel’s and Pearson’s distributions, 
were obtained with the use of both quantile method and 
maximum likewlihood method. The achieved results have 
been presented in Tab. | through 5 and partially in Fig. 1 
through 5. 

— For the determined theoretical sea levels the confidence 
intervals at the confidence levels: Pá = 68,3% 1 Pa = 95,0%, 
were calculated (Tab. 7). 

— Consistency of the assumed theoretical distributions and 
empirical one was checked by means of the Kolmogorow 
test (Tab. 6). Results of the test calcultations do not lead to 
rejection of the assumed hypothesis on consistency of the 
distributions. In this work the Gumbel’s distribution was 
considered the closest to the empirical one. Its advantage is 
lack of lower limitation of variability range of the random 
variable X, among estimated parameters. 

— Among the methods applied in this work the maximum 
likelihood method has appeared the most accurate because 
of lack of necessity of interpolating the distribution 
parameters (in contrast to the case of using the quantile 
method). 


4. 


The obtained results can be considered reliable because of 
the long observation series of sea levels, especially those 
for Gdansk (1887 + 2006), Swinoujscie (1901 + 2006) 
and Kolobrzeg (1867 + 2006). Large (upper) values of the 
observed sea level series are due to storm swellings and 
their impact on sea coast. The occurrence probability of 
high sea levels determined in this work for Swinoujscie, 
Kolobrzeg, Ustka, Wladyslawowo and Gdansk can be used 
in designing the coastal hydro-engineering buildings as well 
as in managing the costal zone and inudation areas during 
storm and flood phenomena. 
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Development of new technologies 
for shipping natural gas by sea 
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ABSTRACT 


In recent years dynamic increase of orders for ships intended for liquified natural gas 
(LNG) shipping has been observed with simultaneous trend of increasing their transport 
capability. This results from the fact that natural gas has become today the third energy 
} source worldwide just next to crude oil and coal. The fast growth of demand for natural 

gas and its limited resources would cause growth of its price, therefore better solutions of 
natural gas transport technology with respect to economy, ecology and safety should be 
searched for. This paper presents various technologies for natural gas transport by sea 


with special attention paid to some alternative methods of transport, namely: CNG and NGH ransport 
technologies in contrast to LNG one. 


Keywords: natural gas; LNG carriers; liquified gas; CNG carriers; compressed gas; 
NGH carriers; gas hydrate; operational cost of gas transport 


INTRODUCTION 


Natural gas industry plays greater and greater worldwide 
role both in economical and political sense. Growth of gas 
consumption is associated with its great popularity as the most 
ecological source of energy. Ecological features of natural 
gas as well as its wide range of applications a.o. to electric 
power production, combustion engine driving and chemical 
production, make the demand for it still growing. The demand 
increases by about 2 % every year. 

The largest natural gas resources are found in the Persian 
Gulf (abt. 41% of worldwide resources) and Russia (abt. 32% of 
worldwide resources), moreover in North Africa and the region 
of the Mexican Gulf and Carribbean Sea, Fig. 1, [2] 
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Fig. 1. Natural gas all over the world [2] 


In Poland natural gas deposits are located mainly in West 
Pomerania, Wielkopolska region and Carpathian piedmont; 
rather small amounts of the gas were also found in sea-bed 
deposits in Polish economical zone of the Baltic Sea. 
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Therefore the increase of number of orders for natural gas 
carriers with simultaneous trend for increasing their transport 
capability reaching today as many as 250000 m°, is observed. 
The boom in production of gas carriers has contributed, due to 
strong competition, to the downward trend in their production 
costs. However it does not change the fact that gas carriers 
belong to very expensive ships because of difficulties in their 
production and operation. 

Natural gas can be shipped by sea in the liquified state 
(LNG), compressed state (CNG), or in the form of hydrate 
(NGH) - the technology being still in the phase of design 
considerations. The CNG and NGH technologies are 
promising, applicable in short-range shipping up to 3000 Mm, 
more profitable, safe and ecological than LNG technology. 
Characteristic features of the natural gas shipped in the three 
forms are presented in Fig. 2. 

A type and design of ship intended for the transporting of 
the raw material as well as equipment of gas loading terminal 
depends on a state in which it will be shipped. 

For many countries, including Poland, gas shipping by sea 
with the use of gas carriers is really the only alternative to be 
independent from gas delivery from Russia and this way to 
greatly improve safety of energy supply to the country. 


CHARACTERISTICS OF LIQUIFIED GAS 
TRANSPORT TECHNOLOGY 


The most popular transport technology of natural gas cargo 
is its liquified form — by means of LNG (liquified natural gas) 
carriers in — 163° C temperature and a little elevated pressure 
of 0.17 MPa (1.7bar). Due to so low temperature LNG must 
be contained in a cryogenic tank. After liquefaction process 


Natural gas shipping by sea 
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Fig. 2. Forms of gas cargo for shipping it by sea and its main properties [the author s elaboration] 


its volume is reduced abt. 600 times against its initial state, 
that constitutes its main advantage for shipping and storing. 
In worldwide trade more than 1/4 of natural gas amount is 
shipped just in this state. 
The LNG is called the concentrated energy as: 
- after regasification, from 1m? of LNG abt. 600 m°? of 
network gas is obtained 
- from 1 t of LNG abt. 1380 m° of network gas is obtained. 


LNG transport chain 


LNG delivery chain begins in natural gas deposits 
from which gas is transported through piping network to 
a liquefaction facility which delivers the gas to an export 
terminal and from this point it is transported, already in the 
liquefied form, by means of special LNG carriers to a import 
terminal and from it, after regasification, it is pumped to gas 
piping network to be delivered to its consumers, Fig. 3. 

According to statistical data from the end of 2008 the world 
fleet of LNG carriers was consisted of 280 ships [10], most of 
which were of over 100 000 m’ capacity. During two last years 
number of built Q-flex and Q-max ships of capacity greater 
than 200 000 m?each, has increased. 

LNG carriers can be divided with respect to their size as 
follows: 

- medmax LNG carriers of the cargo capacity V, = 75 000 m° 
- conventional LNG carriers of the cargo capacity 

V, = 125 000 = 145 000 m° 
- new LNG carriers of the cargo capacity V, = 155 000 = 

+170 000 m’ 


- Q-flex LNG carriers of the cargo capacity V, = 216 000 m° 
- Q-max LNG carriers of the cargo capacity V, = 265 000 m°. 


LNG carriers belong to the group of very specialized ships, 
that results mainly from the applied cargo transport technology. 
The crucial problem is to ensure continuous cooling the cargo 
as well as to avoid its evaporation to the atmosphere. Therefore 
LNG transportation requires application of very effective 
insulation of tanks. 

Among merchant ships LNG carriers are of the highest 
speed, their average service speed varies usually in the range 
of 19 + 21 knots. LNG shipping specificity consisting in 
continuous loss of the cargo during voyage to some extent 
forced voyage period shortening to a minimum. The daily 
Boil-Off-Rate (BOR) amounts on avarage to 0.15 + 0.2 % 
of gas cargo mass and the amount depends mainly on degree 
of insulation effectiveness. The evaporated cargo can be 
discharged to the atmosphere, again liquified or utilized 
onboard as fuel for main propulsion diesel engines, gas turbines 
or boilers delivering steam to turbines. 

The largest LNG carrier (,,Mozah’’) built nowadays, fitted 
with five membane tanks is of 266 000 m’? cargo capacity. 
It was built by Samsung Shipyard in South Korea for 
QGTC company. There are a few other ships like that under 
construction but their size, design or specification have not 
been disclosed so far. 

Attention should be paid to the fact that safety of cargo in 
tanks and ship operation safety should be also ensured along 
with the increased size of the ships. 
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Fig. 3. LNG transport chain with the use of LNG carriers acc. [4] 
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Fig. 4. CNG transport chain by using special gas carriers [the author ss elaboration] 


CHARACTERISTICS OF COMPRESSED 
GAS TRANSPORT TECHNOLOGY 


The CNG is natural gas compressed under high pressure. 
The typical pressure range used to compress the gas is of 
200 + 250 bar. The compresion process is aimed at concentrating 
gas energy in unit volume. The higher pressure the greater 
amount of gas contained in unit volume. After the compression 
about 200/1 + 300/1 units of gas can be achieved. 

Storing and transporting the CNG is carried out in special 
tanks, under high pressure. Capacity of such tanks — usually 
cylindrical vessels — depends on working pressure, temperature 
as well as chemical content of natural gas. The tanks having 
different values of length and diameter, are made of steel or 
light composite materials. 


CNG transport chain with application 
of gas carriers 


CNG delivery chain begins in natural gas deposits from 
which the gas is sent through piping system to compression 
facility and next through loading system (buoy, port) to a ship 
which transports it in the compressed form to reloading terminal 
(buoy, port). In the terminal gas can be reloaded to underground 
stores or directly to gas piping network after passing through 
internal installations of the terminal [9]. Fig. 4 shows the 
schematic diagram of the CNG transport path. 

Concept of gas transportation in the compressed form is 
not entirely new as the first attempts to CNG transporting for 
commercial purposes were made by Columbia Gas Company, 
in the 1960s, but it is one of alternatives for LNG transport, 
especially to servicing local markets and satisfying the needs 
without signing any delivery contracts. 

According to available analyses CNG transport technology 
by sea within 3000 Mm delivery range becomes competitive to 
LNG mode of transport as well as to that by undersea piping 
systems. An advantage of CNG technology is avoidance of 
building an expensive liquefaction facility for natural gas, 
its storing and building tanks to its regasification which are 
often placed close to areas of large density of population, that 
impairs its safety. The new natural gas transport technology 
by using CNG carriers is associated only with building the 
facility for gas compression to over 200 bar pressure in the 
receiving port, that leads to large cost savings within the whole 
gas transport chain. 

During CNG transport the cargo can be reloaded far 
offshore by using loading buoys anchored in a safe distance 
from coast. 

Transport of gas in the compressed form is characterized 
by many advantages, a.o. the following: 

- safety nad ecology — gas loading and reloading can be 
performed far offshore by using loading buoys; in the case 
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of occurrence of a possible leakage CNG does not explode 
but evaporizes and disperses in the atmosphere; when 
combusted it emitts less contaminations to the atmosphere 
than oil or coal, lowering this way emission of noxious 
products such as: CO, CO,, NO, 

- flexibility of gas delivery chains — it means a.o. possibility 
to adjust amount of gas delivery to current demand of the 
market 

- savings in transport and storage costs as well as investment 
cost of the reloading port (abt. 80% of the whole investment 
cost is connected with CNG carrier - Fig. 8). 


Concepts of CNG carriers 
and design solutions of cargo tanks 


In recent years big gas concerns have discovered great 
possibilities in CNG transport technology; many of them have 
elaborated their own CNG transport technologies and designs 
of CNG carriers. The following ones belong a.o. to the most 
important: 

- Coselle (Williams) — gas is transported in turns of steel, 
thick-walled pipes of the diameter D = 0.168 m and of 
1=1700 km in length on average, winded onto the so called 
carousel”, in ambient temperature, under 275 bar pressure. 
The Coselle CNG transport system is modular and intended 
for extension. It means that CNG transport system may 
begin from delivery of gas in small amounts and then be 
developed by adding more and more greater number of the 
gas carriers. Gas cargo capacity of such ships is contained 
in the interval of 1.4 + 56 mln m? within 2000 Mm range 
of delivery [6, 8]. 

According to literature sources the technology is the least 

expensive alternative of CNG delivery which could bring 

30 % cost savings. 


Gas Capacity = 330 mmscf 
Ship Cost = $90 to $120 mil. 


Knutsen OAS — PNG (Pressurized Natural Gas) 
system characteristic of vertical cylindrical tanks (of the 
dimensions: h= 18 +36 m in height, D = 1.04 m in diameter) 
in which gas is transported in ambient temperature under 
250 bar pressure. The gas cargo capacity of such ships varies 
in the interval of 2 + 30 mln m° of compressed gas within 
the 3000 Mm range of delivery [6]. 

Technical details of the concept have not been revealed so 
far hence it is not clear in which way its cargo tanks differ 
from other technologies. 


CETech — a common concept elaborated by three partners: 
Statoil, Teekay - a Canadian firm, and Leif Héegh & Co 
— a Norwegian one. The planned ships will be equipped 
with steel horizontal pipes capable of transporting the 
CNG compressed up to about 200 + 250 bar, in ambient 
temperature. CNG cargo capacity of the largest of them 
will reach 20000 t. 


EnerSea — VOTRANS - the firm has proposed a little 
different method of CNG transport, namely that in which 
gas is compressed to an appropriate pressure (90 + 130 
bar), and next cooled to the temperature in the range of —20 
+ — 40 °C. The so low temperature is applied to achieve 
greater gas density. 

CNG is transported in a set of pipes made of carbon steel, 
placed horizontally, or vertically in the case of smaller 
capacities. The pipes are insulated thermally from the 
surrounding, in relevant tanks. Ship’s cargo compartment is 
consisted of 100 tank modules of the height h = 24+ 36 m. 
Every module contains 24 longitudinal tanks of the diameter 
D = Im each. The tank modules are insulated by cooled 
nitrogen. Cargo capacity of such ships varies in the interval 
of 5.6 + 57 mln m? of CNG, within the 200 + 3000 Mm 
range of delivery [8]. 


Trans Ocean Gas - Cassettes — the firm has proposed 
a unique method of CNG transport consisting in using 
pressure tanks made of glass fibre reinforced plastic (FRP). 
CNG is transported in ambient temperature under 250 bar 
pressure. Cargo- carrying capacity of ships of the concept 
amounts at the most to 28 mln m? of gas. The cargo tank 
is designed in the form of ,,cassettes” inside of which are 
placed 18 vertical pressure cylinders of the dimensions: 
the height h = 12 m and the diameter D = 1.04 m (in three 
groups of 6 cylinders in each). In the cassette of 12 x 12 x 3 
dimensions the natural gas cargo of about 48110 m? in 
volume can be stored. It is approximately equivalent to 
a 10-feet TEU container [1]. Weight of fully loaded cassette 
of standard dimensions can reach about 200 T. An advantage 
of such FRP pressure tanks, as compared with traditional 
steel ones, is their lower weight by abt. 30 % and higher 
resistance to corrosion and cracks. A main disadvantage of 
the technology in question is a very high production cost of 
FRP tank, greater than that of steel tank by 100%. 


E 


AME O 


il 


TransCanada CNG Technologies — in such design 
concept natural gas is transported in long pipes of the 
dimensions: the length | = 12.4 + 36.6 m and diameter 
D = 1.04 m, in a close to ambient temperature, under 
pressure of abt. 210 bar. The pressure tanks made of high 
strength steel are based on CRLP concept. Outer surface of 
the pipe is covered by glass fibre and resin, that contributes 
to a significant drop of its mass.The concept is intended to 
be applied to small ships/barges of cargo capacity of abt. 
2.83 mln m? of gas. 


TransCanada Gas Transport r 
Modules (GTM) CNG concept EA," 


—_ 
-_—_ 


In Tab. 1 examples of the main design parameters of CNG 


carriers of existing concepts are presented. 


CNG ships, resembling bulk carriers, can be built in large 


and medium shipyards as they are not so much sophisticated 
like LNG and LPG carriers are. 
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Tab. 1. Main design parameters of CNG carriers based on different transport technology concepts (acc. the author s elaboration) 


Concepts 
of CNG 
carriers 


Gas volume 
[m] 


EnerSea 


—Votrans ” 6.226.000 


191.12 


Propulsion 
power 
P. [kW] 


Cylinders 
[units] 


20 000 
compres-sed gas 


Knutsen 
OAS 


260 


EnerSea 
-Votrans 
V-800 


267 


Trans 


Canada 2.830.000 


+) 
Coselle 9 339 000 


®, 9 - the ships taken to further considerations. 


108 
,coselles” 


Tab. 2. Values of weight of the tanks and volume of gas shipped in them for two concepts of CNG ships (the author 5 elaboration) 


CNG ship: of EnearSea — Votrans design of Coselle design 


Characteristics 
of cargo tank 


Modules of gas cylinders: 
1 module = 24 cylinders 
Standard dimensions of the tank: 
h=24 + 36 m. D = 1.04 m 
wall thickness t = 0.0335 m - (assumed for 
calculations) 


Diameter: D = 0.168 m 
Wall thickness: t = 0. 0055 m 
Length of pipe turns: 1 = 17702 m 
Total volume: V = 351 m° 
Gas storing temperature: 10°C 
Gas storing pressure: 240 bar 


Weight of steel tank/ 
Weight of pipe turns in one 
„carousel” [t] 


Pa = FID- O-20 HPs 
P= 19.8 + 29.8 


P „= 408 
Weight of „carousel”: 
= abt. 40 
= 448 


carousel 


zb carousel 


p, = 7.801 [t/m ] - density of X-80 HS steel 


Volume of steel tank/ 
Volume of pipe turns in one 
carousel” [m°] 


V 


zb 


V,, = 21.0 +31.2 


4 3 2 
=>nr+nrH 
3 


V„=325 


Weight of compressed gas 
in cargo tank [t] 


P = V 
= 0.22 t/m? 


P spr. gas 


P.,, = 4.6 = 6.9 


* 
zb P spr. ga: 


- mean CNG density 


Gas weight/cylinder weight 
ratio 


n=P,,/P,, = 0.23 


1=P.,/P,, = 0.18 


Weight of gas in all cargo 
tanks on the considered ships 
9% [t] 


P 


gazu 


* 
ny 


P = 4.9 * 1296 = 6350 


P =P 


ws 


*n 


gazu zb 


P = 72 * 108 = 7776 


Weight of all cargo tanks on 
the considered ships ”,*® [t] 


P 


Ww 


= * 
zb Fa ny 


P „= 20.2 * 1296 = 26131 


P 


= * 
wzb P ny 


P „= 408 * 108 = 44064 


Estimated compression ratio 


*) 
_ 6226000 _ 230 


27631 ~ 1 


», 9 - the main design parameters of the ships - as given in Tab. 1. 
V, —natural gas volume (ship load-carrying capacity), 
V — compressed gas volume in cargo tanks (for the ,,Coselle” tank: equal to abt. 325.7 m°; for the gas cylinder: abt. 21.32 m°, at h = 24.4 m) 
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mH) 


~ 35176 7 


270 


T 


Parameters of compressed natural gas 
and gas cargo tanks 


The greatest difficulty in CNG shipping is associated with 
an excessively large tank mass / shipped CNG mass ratio. As 
compared with LNG technology, load-carrying capacity of CNG 
ships is from two to three times lower. Therefore in every CNG 
transport technology above described attempts are undertaken to 
increase volume of transported gas in relation to weight of cargo 
tanks by improving materials used to production of the tanks 
(e.g. application of non-metallic composites, higher strength 
steel) or increasing the mass/volume ratio of gas by searching 
for optimum parameters in pressure-temperature relation. 

To assess volume of a cargo tank and gas to be contained 
in it, appropriate calculations were performed; its results are 
presented in Tab. 2. Out of the presented design solutions 
of cargo tanks the two types have been selected for further 
analysis: the cylindrical steel tank according to the EnerSea- 
Votrans’ concept and the steel tank in the form of thin — walled 
pipes according to the Coselle’s concept. 

As results from the preliminary calculations, despite the 
greater weight of the cargo tanks according to the Coselle 
concept its pressure and temeperature parameters of compressed 
gas result in the greater gas compression ratio, that makes it 
possible to transport the greater amount of gas. 


CHARACTERISTICS OF TRANSPORT 
TECHNOLOGY OF GAS IN THE FORM 
OF HYDRATES 


By using NGH (natural gas hydrate) ships natural gas in the 
form of hydrates can be transported in the conditions of the low 
temeperature T = —20 °C and atmospheric pressure. 


Natt al y 
ji Forming hydrates (NGH) 


Gas deposits 


Gas purification 


Shipping by sea 


Reloading 


Store tank 


Pelletization 
of hydrates 


Gas hydrate is a substance consisted of frozen water 
particles which form cage structures (pellets) in the condition 
of elevated pressure and/or lowered temeperature, in which 
gas particles are kept. 

To form hydrates, presence of water or its vapour 
(in equilibrium state) and hydrate building gases as well 
as an appropriate range of pressure and temperature is 
necessary. 

The conditions for forming hydrates during their production 
and transportation are presented in Fig. 5. 


Hydrate 
Production 


70 


Pressure [bara] 


Hydrate 20 


Transportation 180 m? 


natural gas 


-30 -20 Š -10 0 10 20 
Temperature [°C] 


Fig. 5. Conditions for forming hydrates 
during their production and transportation 


Transport chain of gas in the form of hydrates 


Transport of natural gas by using NGH ships can be split 
into the following phases, Fig. 6: 


Loading 


Store tank 


Gas distribution stations 


Regasification 


Power production plants 


Fig. 6. NGH transport chain (the author 5 elaboration based on [7]) 
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- mining 

- production of NGH in the form of hydrate pulp or 
granulates 

- preparation to transport and loading into ship 

- reloading and regasification — gas recovery 

- delivery to gas provider. 


Design concepts of NGH carriers 


As for many years the transport technology of gas in the 
form of hydrates by using NGH carriers has been still under 
study there are no built ships of such a type. 

Design concepts of NGH carriers are under elaboration 
in a few countries, namely: Japan, Norway, USA and Great 
Britain. Scientific research centres operating in the countries, 
i.e.: Mississippi State University (MSU), Mitsui Engineering 
& Shipbulding Co. Ltd. (MES) or Norwegian University of 
Science and Technology (NTNU) conduct research projects 
on the storing and transporting of natural gas in the form of 
hydrates. 

In 2001 the first design concept of such ship was elaborated, 
Fig. 7, and as results from a long-term feasibility plan, the end 
of the pilotage project is scheduled on the year 2009. 

The NGH ship design concept resembles a traditional 
cargo ship in design, such as tanker or bulk carrier. The ship’s 
hull is characterized by double side and bottom plating and 
subdivision into 12 holds whose tanks are insulated to prevent 
ice forming outside. Loading operation of the hydrate would 
be carried out by using the piping system leading to the tanks. 
During voyage the heat transferred to the cargo through tank 
walls would locally result in a decomposition of the hydrate 
into gas and ice. The so released gas could be utilized as a fuel 
for propelling main engine. 

According to the preliminary design of NGH ship each of 
its 12 tanks was fitted with a mechanical self-reloading device 
similar to those used to reloading operations on bulk carriers. 
The necessity of installing such devices in every cargo tank of 
NGH ship contributes to a very high cost of such ship. 

Transport of gas in the form of hydrates is another 
alternative in relation to LNG transport technology and has 
many advantages, a.o. such as: 

- CNH is easy for storing and safe in transporting 
- low transport requirements as to temperature and pressure 

(T =-20 °C, at atmospheric pressure) 

- low investment and operational costs — inexpensive in 
production and delivery 

- high transported hydrate mass /tank mass ratio 

- lower risk of possible explosion due to cristalline 
structure 


- ecological — not dangerous to the environment 

- possible conversion of existing ships, e.g. oil tankers 

- adjustment of gas deposit resources and quality of 
transported gas to receiver’s needs 

- lower requirements as to thermal insulation of cargo 
tanks. 


However the following is disadvantegous: 
- hydrate contains abt.10% of water 
- small gas volume relative to cubature 
- profitable only for ships of VLCC size 
- necessity of regasification. 


Like in CNG transport technology, only short-range 
shipping up to 3000 Mm is profitable. When properties of the 
gas in LNG and NGH form are taken into account, size of 
the fleet of NGH ships should be almost four times greater as 
compared with that of LNG ships (under the assumption that 
LNG contains 600 m? of natural gas and NGH about150 m° 
of it). Therefore load-carrying capacity of NGH ships must be 
at least four times greater than that of LNG ships at the same 
transported quantity of natural gas, that detrimentally influences 
cost of NGH shipping by sea. 


OPERATIONAL COSTS OF DIFFERENT 
NATURAL GAS TRANSPORT 
TECHNOLOGIES 


The two following factors decide on which technology of 
gas transport by sea would be the most profitable: 
- load-carrying capacity of ship 
- transport range (distance). 


Size of gas carriers as well as their whole fleet and especially 
volumetric capacity of their cargo tanks would first of all 
depend on: 

- direction of import 

- quantity of imported gas 

- land infrastructure, i.e.: gas receiving terminal and gas 
piping network 


whereas their operational costs would be dependent on: 

- location of natural gas import source, i.e. choice of transport 
route 

- fuel cost (depending a.o. on current fuel price, specific fuel 
consumption, ship speed etc) 


275m 
46m 
26m 
14m 


Deadweight: 100,000ton (LNG equivalent :13,000ton) 
Superimposed capacity: 160,000 ni 


Operating speed: 17knot 


Fig. 7. NGH ship design concept [7] 
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- size of ship and kind of its propulsion system (it mainly 
concerns LNG ships) 

- type and number of applied cargo tanks as well as design 
parameters of ship hull 

- personnel cost — freight rate (long-term one) 

- additional expenditures. 


Moreover for calculations of the operational costs the 

following should be assumed: 

- period of active operation of ship 

- period of ship stand-by in ports during each round voyage 
of the ship 

- number of ships angaged on a given route— depending 
on their load-carrying capacity and speed, as well as for 
LNG ships - on a gas cargo quantity evapourized during 
voyage. 


Summing up, operational costs depend on many above 
mentionad factors. Knowledge of all the parameters makes it 
possible to perform optimization investigations to select the best 
variant of gas carrier together with appropriate gas transport 
technology, regarding gas delivery cost minimization. 

Determination of exact operational cost is usually difficult 
and complicated because of lack of reliable data, as well as due 
to commercial secret policy from the side of firms competing 
on the market. 

In Fig. 8, on the basis of literature sources and available 
knowledge, a comparison of component costs of LNG and CNG 
transport technologies is presented, and in Tab. 3 — a comparison 
between costs of LNG and NGH transport technologies. 


Shipping by sea 


Liquefaction 


Reloading 
Costs of application of LNG ships 
Reloading 6% Compression and loading 5% 


Shipping by sea 


Costs of application of CNG ships 


Fig. 8. Comparison of component costs of application 
of CNG and LNG ships [3] 


As results from an economic analysis performed in 2002, the 
costs of the NHG transport chain are lower by 12% in relation 
to those of LNG transport chain. Shipping costs are similar, but 
the option of NGH transport technology is more expensive by 
2% only, like the difference of their regasification costs which 
amounts to 1% only, Tab. 3, acc. [5]. 


CONCLUSIONS 


— Fast growing demand on natural gas and its limited resources 
will cause increasing prices, therefore new, better solutions 
of natural gas transport technology should be searched 
for economical and ecological reasons. Limited possible 
application of piping systems to its transportation forced to 
an extent the development of mass transport technologies 
of natural gas in the liquid or compressed form. 

— Long-range transport of liquified gas is (despite its 
complexity and high cost) the most profitable transport 
method, whereas CNG and NGH transport technologies 
are the best for transportation of medium and small 
amounts of gas over short distances, up to 3000 Mm. 
Therefore economic merits of the two transport alternatives 
should be taken into consideration in planning natural gas 
transport to Poland, especially in the case of its delivery 
from Norway. Expected delivery distance (from Barents 
Sea to Swinoujscie) equal to abt. 2600 km indicates large 
possible savings in investment cost, mainly for transport 
operation. 

— As above mentioned, transport of gas in the form of 
hydrates is still in the research phase hence its practical 
application would be a more distant alternative. For the time 
being, it could be considered a potential, future alternative 
of gas delivery to Poland. 

— Research on development of design and production 
technology of ships intended for natural gas shipping 
are currently carried out in the Department of Ocean 
Engineering and Design of Sea-going Ships, West- 
Pomeranian University of Technology, ZUT in the frame 
of the research project No. 507-09-022-9718/7, titled: 
„METAN — Research on development natural gas shipping 
with special attention paid to the state of transported gas, 
optimum size of ship, type of power plant as well as logistic 
and technolgical problems”. 


NOMENCLATURE 


B — ship breadth [m] 
D — outer diameter of cylinders or pipe turn [m] 
D — maximum ship displacement [t] 
WT — ship deadweight [t] 
h — height of cylindrical tanks [m] 
H — ship hull depth [m] 


Tab. 3. Costs of transport chain based on NGH or LNG technology, acc. [5] 


LNG [%] 
mln $ ([%]) 


Transport chain 


Production 1144 (55%) 


NGH [%] 
min $ ([%]) 


Difference 
mln $ ([%]) 


992 (54%) 152 (13%) 


Shipping by sea with the use of ships 660 (32%) 


628 (34%) 32 (5%) 


Regasification 285 (13%) 


218 (12%) 57 (24%) 


Total 


2089 (100%) 


Note: 
6000 km is equal to 11320000 m’. 


1838 (100%) 251 (12%) 


Calculations were perfromed under the assumption that the yearly amount of gas transported in the delivery range of abt. 
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Psprgas 


r 

t 
Ta 
V 
V, 
V, 
V, 


— length of turn of pipes acc. to Coselle technology [m] 

— ship length b.p. [m] 

— light ship mass [t] 

— number fo cargo cylinders / number of carousels with pipe 
turns [pieces] 

— cylindrical tank diameter [m] 

— cylindrical tank wall thickness [m] 

— maximum ship draught [m] 

— ship speed [m/s] 

— load-carrying capacity of ship [m°] 

— volume of compressed gas in cargo tanks [m°] 

— volume of natural gas [m°] 

— density of compressed natural gas [t/m°] 


ax 
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The Ship Handling Research and Training Centre at Ilawa is owned by the Foundation for Safety of Navigation 
and Environment Protection, which is a joint venture between the Gdynia Maritime University, the Gdansk University 
of Technology and the City of Ilawa. 


Two main fields of activity of the Foundation are: 


> Training on ship handling. Since 1980 more than 2500 ship masters and pilots from 35 countries were trained at 
Iława Centre. The Foundation for Safety of Navigation and Environment Protection, being non-profit organisation 
is reinvesting all spare funds in new facilities and each year to the existing facilities new models and new training 
areas were added. Existing training models each year are also modernised, that's why at present the Centre represents 
a modern facility perfectly capable to perform training on ship handling of shipmasters, pilots and tug masters. 


Research on ship's manoeuvrability. Many experimental and theoretical research programmes covering different 
problems of manoeuvrability (including human effect, harbour and waterway design) are successfully realised at 
the Centre. 

The Foundation possesses ISO 9001 quality certificate. 


Why training on ship handling? 


The safe handling of ships depends on many factors - on ship's manoeuvring characteristics, human factor (operator 
experience and skill, his behaviour in stressed situation, etc.), actual environmental conditions, and degree of water 
area restriction. 


Results of analysis of CRG (collisions, rammings and groundings) casualties show that in one third of all the 
human error is involved, and the same amount of CRG casualties is attributed to the poor controllability of ships. 
Training on ship handling is largely recommended by IMO as one of the most effective method for improving the 
safety at sea. The goal of the above training is to gain theoretical and practical knowledge on ship handling in a wide 
number of different situations met in practice at sea. 


For further information please contact: 
The Foundation for Safety of Navigation and Environment Protection 


Head office: Ship Handling Centre: 
36, Chrzanowskiego Street 14-200 ILAWA-KAMIONKA, POLAND 
80-278 GDANSK, POLAND tel./fax: +48 (0) 89 648 74 90 
tel./fax: +48 (0) 58 341 59 19 e-mail: office@ilawashiphandling.com.pl 

e-mail: office@portilawa.com 
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GDANSK UNIVERSITY OF TECHNOLOGY 


is the oldest and largest scientific and technological academic institution in 
the Pomeranian region. The history of Gdansk University of Technology is 
marked by two basic dates, namely: October 6, 1904 and May 24, 1945. 


The first date is connected with the beginning of the technical education at 
academic level in Gdansk. The second date is connected with establishing of 
Gdansk University of Technology, Polish state academic university. Gdansk 
University of Technology employ 2,500 staff, 1,200 whom are academics. 
The number of students approximates 20,000, most of them studying 


full-time. Their career choices vary from Architecture to Business and (© ¢| d f) S Ki U in} ive rs jy 
of Tacnnology 


Management, from Mathematics and Computer Science to Biotechnology 
and Environmental Engineering, from Applied Chemistry to Geodesics and 
Transport, from Ocean Engineering to Mechanical Engineering and Ship 
Technology, from Civil Engineering to Telecommunication, Electrical and 
Control Engineering. Their life goals, however, are much the same - to 
meet the challenge of the changing world. The educational opportunities 
offered by our faculties are much wider than those of other Polish Technical 
universities, and the scientific research areas include all of 21st Century 
technology. We are one of the best schools in Poland and one of the best 


ESS 
FACul ty) OCEAN huo 
known schools in Europe — one that educates specialists excelling in the = ` pigi 
programming technology and computer methods used in solving complicated E e e 
scientific, engineering, organizational and economic problems. f g | f Id f g 


THE FACULTY OF OCEAN ENGINEERING 
AND SHIP TECHNOLOGY 


The Faculty of Ocean Engineering and Ship Technology (FOEST) as 
the only faculty in Poland since the beginning of 1945 has continuously 
been educating engineers and doctors in the field of Naval Architecture 

and Marine Technology. 


The educational and training activities of FOEST are supported by 
cooperation with Polish and foreign universities, membership in different 
international organizations and associations, as well as participation in 
scientific conferences and symposia. Hosting young scientists and students 

from different countries is also a usual practice in FOEST. 


The activities of Faculty departments are related to: mechanics and 
strength of structures, hydromechanics, manufacturing, materials and system 
quality, power plants, equipment and systems of automatic control, mostly 

in shipbuilding, marine engineering and energetic systems. 


FOEST is a member of such organizations like WEGEMT; The 
Association of Polish Maritime Industries and the co-operation between 
Nordic Maritime Universities and Det Norske Veritas. The intensive teaching 
is complemented and supported by extensive research activities, the core 
of which is performed in close collaboration between FOEST staff and 
industry. We take great care to ensure that the applied research meet both 
the long term and short term needs of Polish maritime industry. FOEST 
collaborates with almost all Polish shipyards. Close links are maintained 
with other research organizations and research institutions supporting the 
Polish maritime industry, such as Ship Design and Research Centre and 
Polish Register of Shipping, where several members of the Faculty are also 

members of the Technical Board. 


The Faculty of Ocean Engineering and Ship Technology is a unique 
academic structure, which possesses numerous highly qualified and 
experienced staff in all above mentioned specific research areas. Moreover, 
the staff is used to effective co-operation and exchange of ideas between 
specialists of different detailed areas. This enables a more integrated and 
comprehensive treatment of research and practical problems encountered 
in such a complicated field of activity as naval architecture, shipbuilding 

and marine engineering. 


The staff of the Faculty has strong international links worldwide, being 
members or cooperating with international organizations like International 
Maritime Organization IMO, International Towing Tank Conference ITTC, 
International Ship and Offshore Structures Congress ISSC, International 
Conference on Practical Design of Ship and other floating Structures PRADS 

just to name a few. 
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